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ABSTRACT 


Future  U.S.  Navy  warships  will  have  DC  electrical  distribution  systems.  These 
distribution  systems  will  include  multiple  layers  of  electronic  power  converters.  When 
coupled  to  high-bandwidth  controllers,  power  converters  behave  as  constant  power  loads 
to  their  distribution  systems.  Constant  power  loads  have  a  negative  non-linear  impedance 
characteristic  that  threatens  system  stability. 

Many  different  single-input  control  schemes  have  been  applied  to  DC  microgrids 
with  constant  power  loads.  In  this  work,  we  propose  a  centralized  select-matrix,  multi-rate 
linear  quadratic  regulator  (LQR-SM)-based  control  scheme  to  provide  a  flexible  and 
adaptable  controller  for  high-order,  multi-input,  and  multi-rate  distribution  systems.  The 
proposed  controller  is  investigated  via  MATLAB  time-domain  simulation.  LQR-SM 
controller  perfonnance  is  compared  to  both  block-cyclic  multi-rate  LQR  and  state- 
feedback  linearization.  LQR-SM  controller  simulations  show  vastly  reduced 
computational  load  and  improved  robustness  compared  to  block-cyclic  LQR  and  reduced 
energy-storage  element  sizing  compared  to  state-feedback  linearization. 

A  genetic  algorithm  design  procedure  is  described  for  the  controller.  The  design 
process  outlines  evaluation  function  development,  choosing  the  initial  generation  of 
candidates,  and  genetic  algorithm  process  flow. 

Finally,  engineering  trade-offs  are  considered.  In  this  work,  we  investigate 
engineering  trade-offs  with  respect  to  computational  load,  regions  of  attraction,  and 
energy-storage  efficiency  when  reduced  fidelity  and  non-adaptive  controller 
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I. 


INTRODUCTION 


Integrated  propulsion  systems  are  the  future  of  naval  surface  combatants.  For 
generations,  U.S.  Navy  ships  have  been  designed  and  built  around  the  concept  of  separate 
power  generation  and  distribution  systems  for  electric  power  and  propulsion  power. 
Demand  for  electrical  power  on  ships  has  been  growing,  fueled  by  the  proliferation  of 
sophisticated  sensors,  communications  equipment,  and  computing  power.  With  the 
additional  development  of  energy  weapons  systems  such  as  lasers  and  naval  railgun, 
demand  for  electrical  power  on  combatant  craft  will  rival  or  exceed  demand  for 
propulsion  power.  If  the  traditional  paradigm  of  separate  electrical  and  propulsion 
systems  is  maintained,  the  physical  size,  weight  and  cost  of  shipboard  engine  rooms  will 
become  prohibitive  [1]. 

Integrated  power  systems  (IPS)  use  a  common  electrical  generation  and 
distribution  system  to  supply  all  loading  on  the  ship.  No  matter  whether  the  loads  are 
typical  hotel  fare,  like  lighting  and  air  conditioning,  energy  weapons,  or  propulsion 
drives,  all  loads  are  powered  from  a  common  distribution  system.  An  IPS  takes 
advantage  of  typical  power  profiles  for  propulsion,  ship’s  service,  and  combat  system 
loads  to  appropriately  size  the  capacity  of  the  generators.  Typically,  full  propulsion 
loading  is  only  used  when  the  ship  is  transiting  between  patrol  areas,  and  low  speeds  are 
used  during  patrol.  Combat  systems  are  not  typically  used  during  transit  periods.  Since 
not  all  loads  will  be  on-line  simultaneously,  these  historical  operating  norms  enable 
electricity  generating  capacity  to  be  smaller  than  total  load  capacity  while  still  meeting 
the  needs  of  the  ship  [2].  IPS  architectures  have  been  installed  on  commercial  cruise 
ships  and  other  vessels,  but  the  first  instance  of  installation  of  an  IPS  on  a  surface 
combatant  is  the  U.S.  Navy  destroyer  USS  Zumwalt  (DDG-1000).  The  USS  Zumwalf  s 
IPS  is  a  4160  V  60  Hz  AC  distribution  system.  While  the  capabilities  of  this  platform  are 
extraordinary,  the  next  generation  of  IPS  warship  is  expected  to  have  a  medium  voltage 
DC  distribution  system.  USS  Zumwalt  has  a  large  number  of  transformers  and  power 
converters  in  order  to  distribute  power  to  meet  various  different  load  needs.  So  far, 
Zumwalf  s  4160  V  60  Hz  AC  power  must  be  converted  to  450  V  60  Hz  AC,  120  V  60  Hz 
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AC,  120  V  400  Hz  AC,  and  low  voltage  DC  power.  The  motivation  for  switching  to  a  DC 
IPS  architecture  is  that  size,  cost,  and  weight  savings  will  be  realized  by  eliminating 
transformers,  reducing  the  size  and  weight  of  cable  runs,  and  minimizing  the  number  and 
type  of  power  conversions.  Furthermore,  DC  distribution  eliminates  the  need  to  match 
frequency  and  phase  when  adding  another  power  source  to  the  distribution  bus.  The  DC 
distribution  also  improves  acoustic  performance  since  generators  and  motors  are  no 
longer  synchronized  to  60  Hz  harmonics.  Additionally,  the  U.S.  Navy  expects  to  find 
performance  gains  by  consolidating  and  centralizing  energy  storage.  Currently,  critical 
devices  are  protected  from  power  outages  by  uninterruptible  power  supplies  (UPS).  Each 
UPS  is  sized  to  meet  worst-case  conditions.  By  moving  to  a  DC  IPS,  energy  storage  can 
be  moved  from  point-of-use  to  points-of-distribution  within  protected  zones.  The 
consolidation  of  UPS  into  fewer  and  larger  energy  storage  devices  should  result  in  lower 
cost,  less  space  consumed,  and  more  efficient  power  conversion  [2], 

A.  MOTIVATION 

A  DC  distribution  system  will  inevitably  have  many  electronics-based  power 
converters.  Electronic  power  converters  are  typically  coupled  with  control  systems  to 
regulate  output  current  and  voltage  or  input  impedance;  thus,  power  converters  do  not 
exhibit  linear  impedance  to  the  distribution  system  in  the  same  way  that  passive  loads  do. 
In  the  worst-case  scenario,  tightly  controlled  power  converters  behave  as  constant  power 
loads  (CPL).  That  is  to  say  that  regardless  of  input  voltage  to  the  power  converter,  the 
power  consumed  by  the  load  is  constant.  CPL  have  negative  non-linear  impedance.  When 
the  net  impedance  of  the  distribution  bus  is  zero  or  negative,  the  distribution  system 
becomes  unstable  [3],  Since  the  future  MVDC  warship  will  very  likely  have  most,  if  not 
all,  of  its  loads  interfaced  and  controlled  through  tightly  regulated  power  converters, 
assuring  stability  of  the  distribution  system  is  a  primary  concern. 

Incorporating  energy  storage  into  the  distribution  system  is  also  a  new  challenge. 
Traditionally,  UPS  and  other  types  of  energy  storage  devices  have  been  used  only  in  “off¬ 
line”  configurations.  The  energy  storage  devices  are  activated  when  a  loss  of  primary 
power  is  detected.  When  primary  power  returns,  the  energy  storage  device  reverts  back  to 
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a  charging/standby  mode.  This  paradigm  for  using  energy  storage  devices  is  inefficient. 
Hybrid  energy  storage  systems  (HESS)  are  energy  storage  devices  that  seek  to  combine 
the  advantages  of  high  power-density  devices  such  as  capacitors  and  high  energy-density 
devices  such  as  batteries  and  flywheels  by  combining  the  devices  with  power  electronics 
and  control  systems  [4],  By  keeping  the  HESS  connected  to  the  power  bus  while  main 
power  is  connected,  the  dynamic  advantages  of  the  HESS  can  be  used  to  regulate  the 
electrical  distribution  system.  In  one  example,  a  HESS  was  used  in  a  hypothetical 
shipboard  power  system  to  cancel  harmonics  introduced  by  wave  motion  on  the 
propulsion  motor  [5].  In  another  example,  HESS  are  a  part  of  energy  management 
schemes  to  allow  total  system  loading  to  exceed  generating  capacity.  In  these  cases, 
whenever  total  electrical  loading  exceeds  generating  capacity,  the  HESS  provide  power 
to  fill  the  deficit  [6],  [7].  While  the  energy  management  schemes  in  [6]  and  [7]  are 
largely  interested  in  using  the  HESS  to  provide  ride-through  power  when  pulsed  loads 
exceed  generating  capacity,  others  have  explored  using  HESS  as  a  means  to  improve 
voltage  regulation  at  all  times.  In  essence,  HESS  become  part  of  the  bus  regulation 
control  strategy. 

B.  OBJECTIVE 

The  primary  goal  of  this  research  is  to  develop  a  control  scheme  to  regulate  the 
bus  voltage  in  a  hypothetical  naval  MVDC  electrical  distribution  system.  The 
hypothetical  naval  MVDC  distribution  system  will  be  derived  from  concepts  published 
by  the  U.S.  Navy’s  Naval  Sea  Systems  Command  (NAVSEA).  While  the  control  scheme 
will  be  applied  to  a  specific  hypothetical  example,  the  scheme  itself  should  be 
generalizable  to  a  wider  class  of  problems.  Furthermore,  we  wish  to  present  a  design 
strategy  for  our  controller  such  that  no  part  of  the  design  process  is  hidden  from  the 
reader.  Any  engineer  should  be  able  to  read  this  document,  follow  the  methods  presented, 
and  obtain  a  satisfactory  result. 

While  previous  authors  have  explored  manifold  methods  to  improve  stability  in 
systems  with  CPL,  we  have  found  their  work  inadequate  for  our  application.  We  have 
discovered  no  sources  within  the  IEEE  Xplore  archive  which  have  explored  a  multi-input 
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approach  to  stabilizing  multi-machine  MVDC  ships  or  microgrids.  We  wish  to  coordinate 
the  actions  of  power  generation  unit  voltages  and  HESS  currents  in  order  to  ensure 
system  stability  under  steady-state  conditions  as  well  as  during  transients  caused  by  large 
step-changes  in  CPL  loading. 

The  secondary  goal  of  our  control  scheme  is  to  minimize  the  volume,  weight,  and 
acquisition  cost  of  the  control  scheme.  As  the  literature  review  shows,  stability  can  be 
assured  through  use  of  various  means.  Many  of  these  means,  such  as  adding  large 
quantities  of  capacitance  or  introducing  passive  resistances,  are  quite  inefficient. 

C.  OUTLINE  OF  DISSERTATION 

The  outline  of  this  dissertation  is  as  follows:  background  infonnation,  such  as  the 
electric  warship  concept,  DC-DC  converters,  average  value  models,  constant  power 
loads,  and  literature  review  of  state-of-the-art  controllers  for  this  application  are 
discussed  in  Chapter  II.  In  Chapter  III,  we  use  the  information  covered  in  Chapter  II  to 
construct  a  hypothetical  circuit  model  for  a  naval  MVDC  IPS  electrical  distribution 
system.  The  proposed  control  method  is  described  in  Chapter  IV.  The  design  approach 
for  the  proposed  controller  is  outlined  in  Chapter  V.  In  Chapter  VI,  we  explore 
performance  trade-offs  incurred  from  different  controller  design  approaches.  Conclusions 
and  future  work  are  presented  in  Chapter  VII. 
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II.  BACKGROUND 


A.  THE  ELECTRIC  WARSHIP 

The  Electric  Warship  [8]  is  a  concept  for  a  naval  surface  warship  for  which  the 
electrical  distribution  system  is  the  central  enabler  of  combat  power.  The  key 
technologies  that  enable  the  Electric  Warship  are  power  electronics,  control  systems,  high 
power-density  electric  motors,  and  energy  weapon  systems. 

Power  electronics  allow  power  to  be  generated  at  one  voltage  and  frequency, 
transmitted  at  another,  and  then  consumed  by  a  load  at  yet  another  voltage  and  frequency. 
By  interspersing  power  electronics  between  the  generation,  transmission,  and 
consumption  stages  in  a  distribution  system,  we  can  optimize  each  different  stage 
according  to  the  needs  of  that  particular  stage.  Advanced  electric  drive  motors  and  motor 
drives  are  key  enablers  since  we  use  these  machines  for  high  power  propulsion  in  our 
electric  ship.  Developing  and  acquiring  high  power-density  machines  and  reliable  drives 
is  absolutely  necessary  for  installation  on  board  a  naval  asset.  Control  systems  tie  the 
elements  of  the  system  together.  Specifically,  the  power  electronics  require  high-speed, 
high  reliability  controllers  to  direct  the  action  of  semiconductor  switches. 

IPS  architecture  integrates  the  propulsion  system  with  the  electrical  distribution 
system.  In  conventional  naval  architectures,  the  electrical  system  and  propulsion  system 
are  independent  systems.  The  propulsion  train  has  prime  movers  connected  to  the 
propulsion  shafts  via  mechanical  transmissions.  In  large  ships,  such  as  aircraft  carriers, 
this  arrangement  has  propulsion  shafts  running  the  entire  1100-ft  length  of  the  ship, 
penetrating  several  watertight  bulkheads.  In  an  IPS,  power  available  from  the  distribution 
system  is  converted  by  a  motor  drive  into  the  fonn  required  by  the  electrical  propulsion 
motors.  The  propulsion  motors  can  be  placed  in  one  compartment,  the  motor  drive  can  be 
placed  in  a  different  compartment,  and  the  power  generating  unit  can  be  in  another. 

The  flexibility  of  an  IPS  has  many  advantages.  First,  prime  mover  speed, 
distribution  system  electrical  frequency,  and  propulsion  motor  speed  are  no  longer 
connected,  improving  acoustic  performance.  Second,  propulsion  and  hotel  loads  can 
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share  power  generating  units.  This  allows  prime  movers  to  be  brought  on  line  or  secured 
to  better  match  generating  capacity  to  loading.  Matching  generating  power  to  loading  can 
be  more  efficient  than  the  standard  arrangement  because  the  gas  turbines  used  for  prime 
movers  are  more  efficient  when  they  are  fully  loaded.  By  amalgamating  loads  onto  fewer 
machines,  greater  fuel  efficiency  can  be  achieved.  Second,  sharing  power  generation  by 
propulsion  loads,  combat  system  loads,  and  hotel  loads  can  reduce  the  number  of  prime 
movers  while  still  maintaining  redundancy.  Reducing  the  number  of  prime  movers 
reduces  acquisition  cost,  maintenance,  and  saves  space  and  weight.  A  third  major 
advantage  is  ship  arrangement  flexibility.  By  eliminating  the  mechanical  connection 
between  prime  movers  and  propulsion  shafting,  we  can  reduce  shaft  lengths  and 
reallocate  compartment  space  for  other  purposes.  Prime  movers  can  be  redistributed  to 
compartments  above  the  waterline  for  better  flooding  protection.  This  allows  exhaust 
stacks  to  be  much  shorter  and  take  up  less  internal  volume  within  the  ship,  as  shown  in 
Figure  1.  In  Figure  1,  we  can  see  that  engineering  spaces  are  designated  in  blue,  freely 
configurable  space  is  shaded  in  brown,  engine  exhaust  paths  are  shaded  in  pink,  while 
propulsion  shafting  is  a  red  line. 


Figure  1.  IPS  Arrangement  Comparison.  Source:  [9]. 
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B.  DC-DC  CONVERTERS 

Power  converters  can  convert  energy  between  alternating  current  and  direct 
current  (AC-DC),  between  one  DC  voltage  level  and  another  (DC-DC),  and  even 
between  one  AC  energy  type  and  another  (AC-AC).  These  electronics  are  largely  based 
on  high-speed,  high-power  rating  semiconductor  switches.  Power  electronics  are  active 
devices  which  require  control  mechanisms  to  regulate  the  switching  action  of  the 
semiconductor  devices. 

The  power  converter  we  discuss  here,  for  illustration  purposes,  is  the  buck 
converter.  The  buck  converter  is  a  DC-DC  converter  that  converts  DC  power  from  a 
higher  voltage  to  a  lower  voltage.  There  are  also  boost  converters  that  convert  energy 
from  lower  voltage  to  higher  voltage,  and  buck-boost  converters  which  are  able  to  both 
buck  and  boost  energy.  Academic  literature  presents  a  wide  variety  of  topologies,  each 
with  its  own  particular  benefits  and  applications.  Rather  than  delve  into  specifics,  we 
review  a  very  basic  buck  converter  since  doing  so  is  instructive  toward  the  development 
of  our  experiment  model.  The  model  of  a  buck  converter  circuit,  also  called  buck 
chopper,  is  shown  in  Figure  2. 


Figure  2.  DC-DC  Buck  Converter.  Source:  [10]. 

When  the  switch  S  is  shut,  the  diode  D  does  not  conduct.  Current  flows  from  the 
voltage  source  E  through  the  inductor  to  the  capacitor  C  and  resistive  load  R.  During  this 
period,  both  the  inductor  and  capacitor  can  build  up  stored  energy.  When  the  switch  S  is 
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open,  no  power  flows  from  the  power  source,  E.  Instead,  energy  for  the  load  comes  from 
the  stored  energy  in  both  the  inductor  and  capacitor.  The  two  equivalent  circuits  for 
switch  “shut”  and  switch  “open”  are  displayed  in  Figure  3.  The  differential  equations  for 
Figure  3a,  the  switch  “shut”  case,  are  given  by 


h~Y  ~L 

Vc  =— - — 

C  CR 


,(2-1) 


and  the  differential  equations  for  Figure  3b,  the  switch  “open”  case,  are  given  by 


Vc 

L 


C 


Vc_ 

CR 


(2.2) 


where  E  is  the  input  source  voltage,  Vc  is  the  capacitor  voltage,  /'/_  is  the  inductor  current, 
and  L,  R,  and  C  are  the  inductor,  resistor,  and  capacitor  values,  respectively.  If  the  period 
when  the  switch  is  “shut”  is  Ton  and  the  period  when  the  switch  is  open  is  T0jf,  the  total 
period  T  is  Ton  +  Toff-  For  a  fixed  period  T,  a  duty  cycle  D  is  defined  as  the  proportion  of 
the  time  the  switch  is  “shut”;  therefore,  D  is  Ton/T.  The  steady  state  analysis  of  [10] 
shows  that  in  the  continuous  conduction  mode  (CCM),  the  output  voltage  Vc  of  the  buck 
converter  is  equal  to  the  input  voltage  multiplied  by  the  duty  cycle  ( Vc  =  DE  ). 


(a) 


(b> 


FIGURE  7.2 


Figure  3.  Buck  Equivalent  Circuits  for  Switch  “Shut,”  (a),  and  “Open,”  (b) 

Source:  [10]. 
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In  order  to  maintain  our  buck  converter  in  CCM,  the  converter  must  have 
sufficient  inductance  for  the  expected  range  of  load.  The  minimum  inductance  the  buck 
converter  needs  to  carry  is  called  the  critical  inductance.  The  formula  for  critical 
inductance  of  a  buck  converter  is 

4,„=  (i-P) 

J switch  Load  n  ti 


where  Lcnt  is  the  critical  inductance,  f switch  is  the  frequency  of  actuating  the  switch  and 
Load  is  the  load  current.  From  Equation  (2.3),  we  see  that  longer  duty  cycles,  higher 
switching  speeds  and  larger  load  currents  will  shrink  the  size  of  our  inductor;  whereas 
larger  loads  and  lower  switching  frequencies  will  do  the  opposite.  We  use  this  critical 
inductance  value  later  when  developing  our  circuit  model. 

The  next  aspect  of  the  buck  converter  we  need  to  explore  is  capacitor  sizing. 
Typically,  switch-mode  power  supplies  will  have  a  specification  for  maximum  voltage 
ripple.  In  order  to  meet  ripple  requirements,  the  filter  capacitor  C  must  be  sufficiently 
large.  The  minimum  capacitor  size  to  meet  ripple  requirements  is 


CT  = 


f\ -£> v 


f 

V  J  switch  J 


100 


8  L  ■  %Vripple 


(2.4) 


where  Cmin  is  the  minimum  capacitor  size,  Vnpp/e  is  the  voltage  ripple  specification,  and  all 
other  parameters  are  as  previously  defined.  In  Equation  (2.4),  we  can  see  that  the 
minimum  capacitor  size  gets  smaller  for  faster  switching  speed,  higher  duty  cycle,  larger 
inductance,  and  larger  ripple  specification. 

C.  AVERAGE  VALUE  MODELS 

In  order  to  examine  the  transient  dynamics  of  the  buck  converter,  we  employ  the 
average  value  model  technique.  Average  value  models  simply  average  the  differential 
equations  of  the  switching  power  converter  over  one  duty  cycle.  Using  this  technique,  we 
may  “average”  the  switch  out  of  the  circuit  and  obtain  a  linear  model  of  the  power 
converter.  In  order  for  our  average  value  model  to  be  valid,  we  must  ensure  that  operation 
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is  limited  to  the  CCM  regime.  In  CCM,  the  current  in  the  inductor  is  always  greater  than 
zero,  thus  inductor  stored  energy  is  never  depleted. 


To  obtain  the  duty  cycle  average  value  model,  we  time-average  the  differential 
Equations  (2.1)  and  (2.2)  over  one  duty  cycle  as  in 


di 
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1  rT  . 

=  —  iL-dr 
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C-AVG 


L  '  C 

C  CR 


(2.5) 


When  we  examine  the  results  of  the  average  value  analysis,  we  see  that  in  steady-state, 
the  capacitor  voltage  is  equal  to  the  duty  cycle  times  the  input  voltage,  and  the  inductor 
current  is  equal  to  the  load  current.  This  result  is  supported  by  the  textbook  analysis  of 
[10].  The  dynamics  of  the  buck  converter  are  exactly  the  same  for  the  dynamics  of  a 
simple  RLC  circuit  but  with  the  input  voltage  modulated  by  the  duty  cycle.  To  complete 
our  model,  we  next  consider  the  current  produced  by  the  voltage  source  E.  By  a  simple 
power  balance,  the  power  produced  by  the  voltage  source  E  (PE  =  ELe  )  must  be  equal  to 

the  power  produced  by  the  average  value  input  voltage  supply  ( PED  =  EDiL  );  therefore, 
Le  =  DiL .  The  full  average  value  circuit  model  is  illustrated  in  Figure  4.  Similar  analysis 

of  elementary  boost  and  buck-boost  converters  yields  similar  average  value  models. 
These  are  presented  for  completeness  in  Figures  5  and  6.  The  main  purpose  of  showing 
the  average  value  models  of  these  various  converters  is  to  demonstrate  that,  regardless  of 
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the  topology,  elementary  switching  power  converters  can  be  reduced  to  simple  second- 
order  models.  Even  though  these  average  value  models  are  simple,  they  preserve  the 
dynamics  of  the  original  system  but  without  the  complexity  of  switch-based  models.  The 
use  of  average  value  models  has  proven  to  be  accurate  and  instructive  in  many 
applications  while  providing  the  benefits  of  using  linear  tools  and  much  simpler  and 
faster  computer  simulations. 


Figure  4.  Buck  Converter  Average  Value  Model 


L 


Figure  5.  Boost  Converter  Average  Value  Model 


L 


Figure  6.  Buck-Boost  Converter  Average  Value  Model 
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D.  CONSTANT  POWER  LOADS 


After  discussing  power  converters,  we  naturally  segue  into  the  concept  of  CPL.  If 
a  power  converter  is  coupled  to  a  control  device,  such  as  a  regulator,  then  the  voltage 
gain  of  the  converter  can  be  adjusted,  via  the  duty  cycle,  to  respond  to  changes  in  loading 
and  input  voltage.  For  example,  in  a  laptop  computer  power  supply,  regardless  of  the 
voltage  coming  from  the  wall  socket,  the  laptop  power  supply  provides  consistent  high- 
quality  DC  power  to  the  laptop  computer.  Any  variation  of  the  input  voltage  or  frequency 
is  counteracted  by  the  regulator  to  ensure  that  the  laptop  receives  the  correct  voltage. 
Since  the  voltage  provided  by  the  converter  is  constant  and  we  assume  that  the  load  is 
also  constant,  the  regulated  converter  consumes  power  at  a  constant  level. 

Naturally,  as  we  have  seen  with  our  average  value  models,  the  power  converter, 
even  with  a  regulator,  has  transient  dynamics;  however,  if  the  transient  dynamics  of  the 
regulated  converter  are  much  faster  than  the  transient  dynamics  of  the  rest  of  the  power 
supply  network  it  is  connected  to,  then  we  may  ignore  the  dynamics  of  the  converter. 
This  is  the  idealized  CPL. 


1.  CPL  Impedance 

An  idealized  CPL  is  a  tightly  regulated  power  converter  whose  dynamics  may  be 
effectively  ignored  [11].  The  idealized  CPL  presents  problems  for  the  design,  analysis, 
and  especially  the  stability  of  the  systems  of  which  they  are  a  part.  This  is  due  to  the  non¬ 
linear  negative  impedance  of  the  CPL.  The  simple  Ohm’s  law  relationship  of  the  CPL  is 
inscribed  as 


p  =  V  I 

1  CPL  r  CPL 1  CPL 


(2.6) 


From  Equation  (2.6),  it  follows  that  the  DC,  or  large-signal,  resistance  of  the  CPL  is 

R, 


V  V  P 

_  '  CPL  _  r  CPL  _  1  CPL 
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1  CPL  M  CPL 


CPL 


(2.7) 
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The  small-signal  resistance  of  the  CPL  is  obtained  by  taking  the  derivative  of  CPL 
voltage  with  respect  to  CPL  current  to  obtain 
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d 


small-signal 


dl, 


-V  = 

r  CPL 


d  R 


CPL 


CPL 


dl  CPL  1 CPL 


1  CPL 


CPL 


V ‘ 


CPL 


CPL 


(2.8) 


From  Equation  (2.8),  we  clearly  see  the  negative  and  non-linear  nature  of  CPL 
impedance.  The  concepts  of  Equation  (2.8)  are  demonstrated  visually  in  the  CPL  I-V 
curve  of  Figure  7.  In  Figure  7,  it  is  apparent  that  CPL  impedance  is  strongly  affected  by 
changes  to  bus  voltage.  As  we  see  in  further  analysis,  the  negative  aspect  of  CPL 
impedance  poses  stability  issues  for  any  DC  distribution  system  that  supplies  CPL.  The 
negative  impedance  is  not  the  most  confounding  aspect  of  CPL  impedance;  instead,  it  is 
the  non-linearity. 

Engineers  have  a  wide  variety  of  tools  available  to  analyze  and  control  linear 
systems.  Time-domain  and  frequency-domain  analysis  tools  are  well-known  and 
understood.  Determining  whether  or  not  a  system  is  stable  is  usually  just  a  matter  of 
determining  the  eigenvalues  of  the  system  matrix  in  the  time  domain  or  checking  gain 
phase  plots  for  sufficient  margin  in  the  frequency  domain.  With  a  non-linear  system, 
stability  analysis  and  control  tools  require  more  advanced  training  and  study. 
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Figure  7.  I-V  Curve  for  96  MW  CPL 


2.  CPL  Stability 

Now  that  we  have  a  model  for  CPL  impedance,  we  can  perform  simple  analyses 
to  gain  insight  on  how  CPL  affect  DC  distribution  system  stability.  For  purposes  of 
illustration,  we  examine  CPL  stability  in  the  context  of  a  system  of  cascaded  buck 
converters.  We  use  the  average  value  model  of  the  buck  converter  developed  previously, 
shown  in  Figure  8.  For  this  relatively  simple  circuit  we  will  first  discuss  stability  in  the 
context  of  a  linearized  system,  and  then  we  will  briefly  touch  on  non-linear  Lyapunov 
stability  analysis. 
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Figure  8.  Average  Value  Model  Buck  Converter  with  CPL 


a.  Linearized  System  Analysis 

For  the  system  of  Figure  8,  we  first  derive  the  differential  equations  from 
Kirchoff  s  voltage  law  and  Kirchoff  s  current  law  to  obtain 


lG 


-R  .  1  1 

L  lC  i  Bus  +  i 


1 


V  =  —i  - 

r  Bus  ^  lG 


1 


c  CRcpl 


-vD, 


(2.9) 


where  RCPL  is  the  linearized  CPL  impedance  and  all  other  parameters  are  as  previously 
defined.  Next,  we  use  the  characteristic  polynomial  to  establish  stability  conditions.  We 
obtain  the  characteristic  polynomial  by  first  re-writing  the  differential  equations  in  state- 
space  form,  then  taking  the  determinate  as  in 
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CPL  J 


The  purpose  of  this  exercise  is  to  gain  insight  into  how  changing  certain  circuit 
parameters  will  affect  stability.  The  first  thing  we  must  recall  is  that  Rcpl  is  negative.  In 
order  for  the  system  described  by  Equation  (2.10)  to  be  stable,  all  of  the  tenns  must  have 
the  same  sign.  When  examining  the  second  tenn  of  the  characteristic  polynomial  in 
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Equation  (2.10),  we  see  that  increasing  R,  reducing  L,  increasing  C,  and  increasing  the 
magnitude  of  Rcpl  all  improve  stability.  In  regard  to  these  methods,  increasing  R  is 
undesirable,  since  this  reduces  the  steady-state  power  efficiency  of  the  system.  The 
inductance  L  is  usually  reduced  to  minimum  values  anyway  since  inductors  have  parasitic 
resistance  and  are  generally  large  components.  As  we  saw  in  part  B  of  this  chapter,  there 
is  a  limit  to  how  small  L  can  be;  it  must  be  equal  to  or  larger  than  the  critical  inductance. 
Increasing  the  magnitude  of  Rcpl  can  be  achieved  by  raising  system  voltage  or  by 
reducing  load  power.  This  leaves  increasing  bus  capacitance  as  the  free  design  variable 
for  improving  stability  margins.  Depending  on  the  voltage  and  power  of  the  system,  we 
find  that  bus  capacitance  levels  can  be  quite  large,  with  values  in  the  tens  of  mF. 

Despite  all  of  this  analysis,  we  must  recall  that  this  is  only  a  linearization  of  the 
system.  This  analysis  is  only  valid  for  small  deviations  about  the  nominal  Vbus  voltage. 
Any  large  deviation,  such  as  one  that  might  be  induced  by  a  large  load  transient  or  a  fault, 
could  violate  our  assumptions  and  result  in  operation  within  an  unstable  region. 

b.  Lyapunov  Analysis 

Lyapunov  analysis  techniques  are  required  to  assess  large-signal  regions  of 
attraction  (ROA).  We  follow  the  stability  analysis  methods  of  [12]  to  obtain  the  necessary 
conditions  for  the  circuit  of  Figure  8.  First,  we  must  transform  the  first  order  differential 
equations  of  Equation  (2.9)  by  a  change  of  variables  where  (Vo,  Io)  represents  the  desired 
equilibrium  operating  point.  We  define  the  new  state  variables  as  Z;  and  Z?  according  to 


Z  =  V  -V 

^1  V  Bus  V  0 
■^2  =  lG  ~  A) 


(2.11) 


We  substitute  the  new  state  variables  from  Equation  (2.11)  into  Equation  (2.9)  to  obtain 


f=-f(Zi+/„4(Z,+ro)  +  f 
^  =  i(z2+i0) — — — - 

dt  Cv  2  0J  C(Z,+V0) 


(2.12) 
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If  we  assume  E  =  RI0  +  V0 ,  Equation  (2. 12)  can  be  rewritten  as 


Zi+ - 

1  LC 


RC  — 


LP 


(Zi+v0y  j 


z,  + 


1  LC 


z  _(Ri0z ,) 


+  K0)y 


=  o 


(2.13) 


Since  Equation  (2.13)  is  a  second-order  equation  of  the  form  of 

Z\+a-f(Zl)Zl+b-g(Zl)  =  0 


(2.14) 


a  Lyapunov  function,  \j/,  can  be  found  according  to 


^  =  b-Jg  (P)dP  +  ^~ 
0  ^ 

¥  =~a‘f(Z1)Z^ 


(2.15) 


The  resulting  Lyapunov  function  and  Lyapunov  derivative  are 


¥ 


1 


CbusLgl 


'^-RgU^-RgU^log 


rz1+v^ 

v  K  J 


i  ,  70z, 


ICbus2  j  I  2  ZI  +  K0J  j(216) 


and 


v 


LglCbus  J 


f 

RglCbus 

V 


LglP 


z, 


(2.17) 


respectively.  Lor  an  initial  condition  (v,/)  to  be  within  the  ROA  centered  about  the 
equilibrium  point  ( Vo,  Io ),  the  Lyapunov  function,  Equation  (2.16),  must  be  positive 
definite,  and  its  derivative,  Equation  (2.17),  must  be  negative  definite.  We  substitute 
Equation  (2.11)  back  into  Equations  (2.16)  and  (2.17)  and  obtain  the  inequalities 
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La  i 


2  Cbus 


{i-Io  +  Io*  ^}2  >  {-  \  (v  -  Vo)2  +  Rgl  *l0[v-V0-  V0  In  (^)]}(2.18) 


and 


v-sqrt{L^^)' 


(2.19) 


The  inequalities  of  Equations  (2.18)  and  (2.19)  detennine  the  ROA  for  the  system 
and  chosen  input.  Initial  conditions  within  the  ROA  eventually  settle  back  to  equilibrium 
conditions.  Those  initial  conditions  that  do  not  satisfy  Equations  (2.18)  and  (2.19)  result 
in  oscillatory  or  unstable  system  dynamics.  An  example  ROA  is  presented  in  Figure  9  for 
the  system  of  Figure  8.  Capacitance  is  varied  to  show  the  effects  bus  capacitance  has 
on  the  ROA.  The  region  displayed  is  limited  to  a  region  from  zero  volts  to  double 
the  steady-state  voltage  and  from  zero  amperes  to  double  the  steady-state  value.  From 
Figure  9,  we  see  that  the  ROA  is  an  ellipsoid  centered  about  the  equilibrium  point  whose 
major  axis  is  oriented  somewhat  parallel  to  the  current  axis  and  whose  minor  axis  is 
parallel  to  the  voltage  axis.  By  inspection  of  this  ROA,  we  can  conclude  that  the  system 
is  far  more  sensitive  to  changes  in  voltage  than  to  changes  in  current.  Furthermore,  small 
changes  in  capacitance  can  have  enormous  consequences  for  the  size  of  the  ROA.  While 
this  analysis  is  valid  for  the  circuit  of  Figure  8,  conducting  non-linear  analysis  for  high- 
order  systems  can  be  quite  complicated.  Finding  appropriate  Fyapunov  functions  is  more 
art  than  science.  While  some  Fyapunov  functions  are  known  to  be  appropriate  for  certain 
systems,  any  change  of  the  dynamics  of  that  system  can  render  the  Fyapunov  function 
invalid. 
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Figure  9.  ROA  for  a  12  kV  Second  Order  System  with  96-MW  CPL 

E.  RELATED  WORK 

A  great  number  of  authors  have  applied  control  techniques  to  expand  the  ROA  in 
DC  power  systems  with  CPL.  The  controls  applied  to  this  problem  include  a  cornucopia 
of  techniques.  For  ease  of  understanding,  we  summarize  the  control  techniques  applied  to 
this  problem  into  two  main  groups:  linear  and  non-linear. 

1.  Linear  Methods 

Linear  methods  are  those  methods  developed  for  control  of  linear  systems. 

a.  Passive  Damping 

Passive  damping  techniques  involve  linearizing  the  CPL  about  an  expected 
equilibrium  point.  From  this  equilibrium  point,  Equation  (2.10)  can  be  used  to  select 
appropriate  component  values  for  the  desired  performance.  This  method  can  be  reduced 
to  minimizing  inductance,  using  large  capacitors,  placing  large  resistances  in  parallel 
with  the  CPL,  or  placing  a  damper  resistance  in  line  with  the  capacitor  as  suggested  by 
Hodge,  Flower,  and  Macalinden  [11]. 
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b.  Active  Damping 

Active  methods  improve  dynamics  by  including  some  variant  of  a  proportional, 
integral,  and  derivative  (PID)  controller  to  provide  the  desired  dynamics.  These  methods 
are  sometimes  dubbed  “virtual  impedance”  methods.  Magne  et  ah,  in  [13],  apply  the 
control  signal  to  the  CPL  regulator  to  introduce  a  “virtual  capacitor.”  A  passivity  based 
criteria  is  used  to  develop  the  PID  gains  in  [14].  Zadeh  et  al.  use  standard  PID  control 
applied  to  both  a  current  control  loop  and  a  voltage  control  loop  in  [15]  and  [16]  to 
stabilize  system  voltage.  Grillo  et  al.  [17]  apply  proportional  current  and  voltage  gains  for 
his  controller. 

The  approach  of  Carmeli  et  al.  in  [18]  is  interesting  in  that  rather  than  controlling 
the  supply  voltage  or  the  CPL  controller  response,  they  chose  to  place  an  energy  storage 
device  (ESD)  in  parallel  with  the  CPL.  They  then  use  a  PI  controller  on  the  ESD  to 
source  or  sink  current  to  enhance  the  stability  of  the  system.  In  the  conclusions,  [18] 
notes  that  compared  to  a  bulk  capacitor,  the  ESD  provides  a  better  dynamic  response  as 
well  as  providing  the  opportunity  to  store  the  energy  at  lower  voltage  than  the  main  bus 
voltage.  This  is  especially  important  for  MVDC  as  capacitors  rated  in  the  10  kV-20  kV 
range  are  heavier,  bulkier,  costlier,  and  more  failure  prone  than  equivalent  capacitors 
rated  at  lower  voltages. 

2.  Non-linear  Methods 

a.  Adaptive  State-Feedback 

Arcidiacono  et  al.  [19]  demonstrate  use  of  an  adaptive  proportional  state-feedback 
controller  to  ensure  stable  operation  throughout  the  full  range  of  operation  for  a  second- 
order  system.  The  proposed  controller  adaptively  updates  proportional  voltage  and 
current  feedback  gain  to  control  the  duty  cycle  of  a  supply  DC-DC  converter  in  order  to 
provide  stable  and  consistent  regulation. 

b.  Feedback  Linearization 

Feedback  linearization  is  a  popular  technique  that  involves  a  two-part  control 
process.  In  the  first  part  of  the  control  process,  whenever  possible  the  state  of  the  system 
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is  properly  redefined  so  that  the  resulting  system  is  linear.  The  second  part  of  the  process 
involves  using  traditional  linear  pole-placement  techniques,  such  as  PID  controllers,  to 
obtain  the  desired  system  behavior.  Ciezki  and  Ashton  [20]  applied  this  technique  to  a 
buck  converter  with  an  attached  CPL  in  the  late  1990s.  Sulligoi  et  al.  developed 
feedback-linearization  controller  implementations  in  [21]  and  [22]  for  multi-machine 
multi-CPL  systems.  Their  approach  reduced  a  multi-machine  problem  to  a  second-order 
single-input  representation.  It  was  further  demonstrated  in  [23]  that  feedback 
linearization  can  provide  stability  under  controller  saturation  conditions.  In  [24],  active 
damping  control  is  effectively  used  for  the  generating  DC-DC  converter  controller 
while  feedback  linearization  is  used  for  main  bus-connected  buck  converters  serving 
CPL.  This  method  was  proved  to  be  far  superior  to  active  damping  when  disturbances 
were  large  [25]. 

c.  Linear  Quadratic  Gaussian 

Linear  Quadratic  Gaussian  (LQG)  control  is  a  technique  where  there  is 
uncertainty  in  the  state  of  the  system  due  to  white  Gaussian  noise  and  modeling  errors. 
The  state  variables  are  estimated  using  a  Kalman  filter  (EKF)  and  a  linear  quadratic 
regulator  (LQR)  is  then  used  to  provide  state-feedback  gains  that  minimize  a  quadratic 
cost  function.  Zhu,  Liu,  and  Cupelli  used  this  technique  in  [26]  to  provide  decentralized 
control  of  a  multi-machine  multi-load  hypothetical  MVDC  shipboard  system  with  CPL. 
By  using  the  EKF,  the  authors  were  able  to  represent  each  generating  DC-DC  converter 
as  having  a  single  CPL,  thus  reducing  the  order  of  the  problem  to  2nd  order.  This 
technique  does  not  require  the  significant  number  of  simplifying  assumptions  needed  in 
[22]  to  reduce  a  complex  system  to  2nd  order.  Load  sharing  between  machines  is 
enforced  by  proportioning  the  estimated  load  calculated  by  the  EKF.  In  this  way, 
robustness  on  loss  of  generating  power  is  accomplished  by  sensing  the  loss  and  re¬ 
apportioning  load  to  the  remaining  generators.  Performance  of  LQG  was  validated  in  a 
hardware-in-the-loop  simulation  [27]. 
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Neither  [26]  nor  [27]  describe  the  construction  of  the  EKF  or  the  Q  and  R 
matrices  which  describe  model  and  observation  uncertainties  required  for  LQG 
implementation.  This  omission  confounds  this  author’s  ability  to  reproduce  their  work. 

d.  Lyapunov  Backstepping 

Backstepping  is  a  technique  where  a  non-linear  control  law  is  developed  using  a 
chosen  Lyapunov  function.  Since  CPL  power  levels  are  variable  in  real-life  applications, 
these  values  must  be  estimated.  Adaptive  backstepping  introduces  estimator  variables  as 
well  as  other  design  variables  necessary  to  ensure  Lyapunov  stability  conditions.  Like 
[22]  and  [26],  Cupelli  et  al.  [28]  again  reduce  the  system  to  second  order  so  that  a  known 
Lyapunov  function  can  be  used.  Backstepping  produces  controllers  which  outperform 
LQG  in  terms  of  overshoot,  undershoot,  and  settling  time  response  as  shown  in  both  [28] 
and  [29].  Unfortunately,  the  claimed  performances  depend  upon  parameter  choices  which 
are  not  well-defined.  As  mentioned  previously,  this  method  is  only  applicable  in  cases 
where  the  system  Lyapunov  function  is  known. 

e.  Synergetic  Control 

Synergetic  control  is  a  method  where,  after  deriving  the  state  space  model,  macro 
variables  are  defined  from  which  control  laws  are  calculated.  This  method  was  applied  by 
Kondratiev  and  Dougal  [30],  [31].  Cupelli  et  al.  [32]  found  synergetic  control 
performance  superior  to  feedback  linearization;  however,  synergetic  control  is  much 
more  difficult  to  grasp  than  previously  mentioned  methods  and  has  not  gained  popularity. 

F.  SUMMARY 

In  this  chapter,  we  explored  the  concept  of  the  Electric  Warship.  To  make  the 
Electric  Warship  a  possibility,  the  U.S.  Navy  is  pursuing  the  development  of  MVDC 
distribution  systems.  These  distribution  systems  rely  on  the  perfonnance  of  a  wide 
variety  of  electronic  power  converters. 

Switch-based  power  electronic  devices  are  inherently  non-linear;  however,  in 
systems  where  the  device  is  switched  according  to  a  duty  cycle  over  a  defined  period, 

linear  average  value  models  of  the  power  converters  may  be  developed.  Linear  average- 
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value  models  simplify  the  analysis  and  simulation  of  power-electronic  converters  and 
provide  a  familiar  entry  point  for  developing  control  strategies. 

The  use  of  power  converters  that  are  tightly  regulated  by  high  bandwidth 
controllers  result  in  CPL.  CPL  present  a  non-linear  negative  impedance  characteristic, 
which  complicates  stability  analysis  and  control  strategies.  While  a  system  may  be  stable 
at  an  equilibrium  point,  the  system  itself  is  not  globally  stable.  Unlike  a  linear  system  that 
is  either  stable  or  unstable,  non-linear  systems  are  only  stable  within  an  ROA  about  the 
equilibrium  point.  Further  complicating  matters  is  that  Lyapunov  analysis  for  high-order 
systems  can  be  very  difficult  or  impossible,  since  finding  appropriate  Lyapunov  functions 
is  more  art  than  science. 

Various  control  strategies  have  been  employed  to  improve  ROAs  for  systems  with 
CPL.  The  earliest  methods  explored  were  passive  damping  methods,  which  add 
additional  resistance  and  capacitance  to  the  system.  Passive  damping  methods  are  simple 
and  reliable,  but  they  introduce  weight,  cost,  and  inefficiencies.  Linear  control  methods 
employ  traditional  PID  controllers  to  expand  the  ROA.  Linear  control  methods  are  more 
complicated  than  simple  passive  damping  but  utilize  the  massive  body  of  knowledge 
available  to  engineers  for  control  of  linear  systems.  Unfortunately,  the  system  is  non¬ 
linear,  so  global  stability  is  not  guaranteed. 

Non-linear  control  methods  have  also  been  employed.  Some  non-linear  methods, 
such  as  adaptive  linearization  and  state-feedback  linearization,  are  fairly  straightforward 
and  user-friendly.  This  is  due  to  their  close  similarity  to  linear  methods.  Strict  non-linear 
methods,  such  as  adaptive  backstepping  and  synergetic  controls,  are  available  but  have 
found  little  popularity  due  to  complexity  and  ambiguity  in  developing  control  laws. 

With  the  exception  of  the  obtuse  synergetic  control,  all  of  the  control  schemes 
surveyed  are  single-input  schemes.  Control  schemes  in  [15],  [  1 9]— [2 1  ]  approach  multi¬ 
machine  systems  by  either  simplifying  multiple  co-located  input  devices  into  a  single 
device  or  by  decomposing  the  system  into  multiple  independent  systems  with  a  single¬ 
input.  The  approach  of  simplifying  co-located  input  devices  into  a  single  input  device  is 
inadequate  for  systems  with  distributed  input  devices.  Neither  does  that  approach  allow 
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input  devices  that  are  significantly  dissimilar.  The  second  approach,  used  by  the  LQG 
controllers,  decomposes  the  multi-input  system  into  several  independent  single-input 
mini-systems.  This  approach  works  well  for  distributed  systems  and  where  dissimilar 
input  devices  are  used.  A  disadvantage  of  that  approach  is  that  controllers  act 
independently  to  regulate  system  performance.  The  decentralized  and  uncoordinated 
nature  of  this  approach  introduces  the  possibility  that  the  separate  controllers  may 
contradict  one  another,  resulting  in  suboptimal  regulation. 
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III.  NAVSEA  CIRCUIT  MODEL 


Before  we  explore  an  appropriate  controller  for  a  shipboard  electrical  distribution 
system,  we  first  need  a  model  for  it.  The  starting  point  for  this  analysis  is  a  proposed 
architecture  for  a  next  generation  Electric  Warship  from  [33].  In  this  chapter,  we  walk 
through  the  development  of  an  experimental  circuit  model.  In  the  first  section,  we 
describe  the  key  elements  of  the  distribution  system.  Next,  the  NAVSEA  proposed 
distribution  system  is  simplified  into  an  abbreviated  system.  Finally,  an  equivalent  circuit 
structure  is  developed. 

A.  NAVSEA  PROPOSED  DISTRIBUTION  SYSTEM 

A  proposed  naval  electrical  distribution  system,  described  in  [8]  and  illustrated  in 
Figure  10  provided  by  Dr.  Norbert  Doerry  of  NAVSEA,  has  a  MVDC  main  distribution 
system.  While  several  voltage  ratings  are  possible,  12.0  kV  appears  to  be  the  most  likely 
candidate. 

The  concept  includes  four  power  generation  modules  (PGMs)  which  can  be  either 
gas  turbine  or  diesel  engine  prime  movers  driving  split-winding  multi-phase  AC 
generators.  The  output  of  the  AC  generators  are  rectified  and  used  as  the  supply  voltage 
for  a  DC-DC  converter.  In  the  proposed  NAVSEA  system,  there  are  two  large  PGMs, 
called  Main  Turbine  Generators  (MTG),  and  two  smaller  PGMs,  called  Auxiliary 
Turbine  Generators  (ATG).  The  four  PGMs  are  connected  to  Port  and  Starboard 
distribution  buses. 

The  distribution  buses  span  several  zones,  each  zone  representing  a  watertight 
boundary  or  functional  grouping  within  the  ship.  The  zones  contain  switchboards, 
intennediate  DC-DC  power  converters,  ESD,  and  point-of-load  DC-DC  and  DC-AC 
power  converters.  The  point-of-load  power  converters  transform  the  intermediate  DC 
power  from  the  intermediate  DC-DC  converters  to  low-voltage  DC  and  the  various  forms 
of  AC  power  used  within  a  ship:  120  V/240  V/450  V  60  Hz  and  120  V  400  Hz. 

The  load  types  can  be  divided  into  three  main  groups.  The  first  group  is  referred 

to  as  Ship’s  Service  loads.  These  are  loads  such  as  lighting,  air  conditioning,  wall 
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sockets,  galley  equipment,  and  other  various  small  loads.  Ship’s  Service  loading  is 
generally  considered  constant. 
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Unpublished  diagram  provided  by  Dr.  Norbert  Doerry  of  NAVSEA 


Figure  10.  NAVSEA  Proposed  Electrical  Distribution  System. 
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The  next  type  of  load  is  referred  to  as  Combat  System  loads.  These  loads  are 
associated  with  energy  weapon  systems.  Whether  the  load  is  a  high-power  radar,  laser 
weapon  system,  or  naval  railgun,  these  loads  are  characterized  as  pulsed  or  stochastic 
type  loads.  Combat  System  loads  can  experience  large  step-changes  in  power 
consumption.  A  free-electron  laser,  for  example,  can  consume  power  on  the  order  of 
10  MW  during  discharge  but  consumes  much  less  power  during  standby  conditions. 
Likewise,  a  naval  electromagnetic  railgun  (EMRG)  consumes  a  relatively  small  amount 
of  power  in  standby  compared  to  the  period  when  it  is  charging  its  energy  storage  system. 
An  EMRG  operating  at  maximum  cyclic  rate  of  fire  could  require  25  MW  or  more  [34], 

The  final  class  of  load  is  propulsion  motors.  Propulsion  loads  have  historically 
been  the  largest  loads  on  a  ship.  The  U.S.  Navy  Zumwalt  (DDG-1000)  class  destroyer  is 
a  15000+  ton  warship  with  two  advanced  induction  propulsion  motors  rated  at  35.0  MW 
each  for  a  total  of  70.0  MW  of  propulsion  capacity  [35]. 

B.  SIMPLIFIED  DISTRIBUTION  SYSTEM 

To  aid  in  producing  an  experimental  circuit  model,  we  reduce  the  diagram  shown 
in  Figure  10  into  a  simplified  model.  The  simplified  distribution  system  should  retain  the 
key  features  shown  in  Figure  10.  The  first  key  feature  that  we  notice  is  the  multiple 
PGMs.  These  PGMs  come  in  two  capacities:  MTG  and  ATG. 

The  second  feature  we  extract  is  multiple  power  conversion  layers.  The  first  level 
of  power  conversion  is  from  the  PGM  to  the  main  DC  buses.  The  second  level  is  from  the 
main  DC  buses  to  the  zonal  DC  buses  via  the  intermediate  DC-DC  converters.  Third 
level  and  greater  power  conversion  layers  occur  via  point-of-load  converters  within  the 
zones.  Since  we  wish  to  focus  our  attention  on  stability  and  control  of  the  MVDC  main 
buses,  we  simplify  all  third-level  and  beyond  power  converters  by  representing  them  as 
ideal  constant  power  loads.  This  simplification  is  reasonable  since  these  converters 
operate  at  lower  voltages  and  can  be  cost-effectively  switched  at  very  high  speeds  (i.e., 
/switch  >  20  kHz).  If  we  make  the  additional  assumption  that  the  point-of-load  converters 
are  tightly  regulated,  we  have  the  necessary  conditions  to  ascribe  CPL  behavior  to  the 
third-level  and  greater  power  converters. 
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The  third  feature  we  extract  from  Figure  10  is  the  zonal  architecture.  We  construct 
our  model  with  three  zones  based  on  the  three  types  of  loads  present.  Zone  #1  contains 
Ship’s  Service  loads,  which  make  up  a  relatively  small  fraction  of  overall  loading  (<20%) 
and  have  relatively  small  load  variation.  Zone  #2  contains  Combat  System  loads,  which 
are  pulse-type  loads,  capable  of  huge  jumps  in  power  in  a  matter  of  microseconds  (up  to 
25-MW  pulses).  Zone  #3  holds  the  propulsion  loads,  which  are  the  largest  loads  on  the 
ship.  Propulsion  loading  can  be  either  static  or  highly  variable  depending  on  the  time- 
scale.  Maneuvers  from  all-back  full  to  ahead  flank,  as  well  as  sea-swell  periods,  are  on 
the  order  of  seconds.  In  order  to  test  worst-case  conditions,  we  assume  that  all  loads  step 
instantaneously. 

The  final  feature  we  extract  is  the  placement  of  ESD.  When  we  look  carefully  at 
Figure  10,  we  notice  that  every  intermediate  DC-DC  converter  has  a  HESS  unit 
associated  with  it;  therefore,  the  Ship’s  Service  and  Combat  System  zones  have  an  ESD 
connected  to  their  intermediate  buses.  Propulsion,  however,  does  not  have  an  associated 
HESS,  since  the  ship’s  inertia  is  sufficient  for  stored  energy. 

The  culmination  of  our  simplifications  is  illustrated  in  the  block  diagram  of 
Figure  11.  Figure  11  has  all  of  the  key  features.  The  four  PGMs  are  two  different  sizes. 
The  two  larger  PGMs  are  40.0-MW  capacity  while  the  two  smaller  PGMs  are  10.0-MW 
capacity  for  a  total  generating  capacity  of  100.0  MW.  The  MVDC  bus  is  nominally  set  at 
12.0  kV  DC.  There  are  three  load  zones.  The  Ship’s  Service  intennediate  DC-DC 
converter  transfers  the  12.0-kV  MVDC  bus  power  to  a  1.0-kV  DC  Ship’s  Service  bus. 
The  Ship’s  Service  zone  is  a  20.0-MW  maximum  load  CPL  with  a  bus-connected  ESD. 
The  Combat  System  intermediate  DC-DC  converter  transfers  MVDC  bus  power  to  a 
6.0-kV  bus.  The  Combat  System  zone  has  a  30.0-MW  maximum  capacity  CPL  with  a 
bus-connected  ESD.  The  Propulsion  intennediate  converter  transfers  12.0-kV  MVDC 
power  to  10.0  kV  DC  for  an  80.0-MW  maximum  capacity  CPL.  It  is  important  to  note, 
that  although  maximum  possible  loading  exceeds  maximum  available  generating  power, 
exercising  the  plant  through  overload  conditions  is  beyond  the  scope  of  this  document. 

A  harmonic  filtering  element  is  included  at  the  input  to  the  intermediate  DC-DC 

converters  since  such  filters  are  often  necessary  in  practical  applications. 
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Figure  11.  Simplified  Distribution  System  Model 


C.  EXPERIMENTAL  CIRCUIT  MODEL 

Transforming  the  distribution  system  model  of  Figure  11  into  a  circuit  model  is 
our  next  task.  We  systematically  move  from  the  PGMs  to  the  loads  in  assigning  a  circuit 
model  to  each  block  in  the  diagram. 

PGMs  are  the  first  blocks  we  tackle.  The  PGMs  are  composed  of  a  prime  mover, 
multi-phase  AC  generator,  rectifier,  and  DC-DC  converter.  In  order  to  make  our  problem 
tractable  and  relatively  easy  to  simulate,  we  assume  that  the  rectified  AC  from  the 
generator  is  an  ideal  DC  source.  While  not  strictly  true,  it  allows  us  to  ignore  motor  and 
generator  rotational  dynamics.  This  is  a  desirable  assumption,  since  developing  a  more 
realistic  model  for  the  generator  would  be  quite  time-consuming  and  ultimately  have  little 
to  no  bearing  on  the  results  of  this  research.  This  assumption  also  allows  us  to  assign  the 
PGM  the  form  of  the  DC-DC  converters  in  Chapter  II,  Figures  4-7.  We  can  further 
simplify  the  circuit  model  for  the  PGM  if  we  simply  ignore  the  left-hand  side  of  Figures 
4-7.  This  reduces  the  PGM  circuit  model  into  a  controlled  voltage  source,  a  series  RL 
impedance,  and  a  parallel  filter  capacitor.  This  way,  we  do  not  need  to  know  what  type  of 
DC-DC  converter  is  being  used  for  the  PGM  DC-DC  converter.  Since  the  PGMs  are 
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operating  at  a  fairly  high  voltage  (12.0  kV),  we  assume  a  DC-DC  converter  switching 
speed  of  1.0  kHz.  This  speed  is  easily  possible  with  existing  semiconductor  technology. 

Next,  we  examine  the  bus  cabling.  Using  the  short-wire  impedance  model  for  the 
MVDC  bus,  we  represent  MVDC  bus  cabling  as  a  series  RL  impedance  with  a  parallel 
capacitance.  This  model  has  basis  in  other  work,  as  seen  in  [27].  For  the  sake  of 
complicating  the  model  such  that  the  simplifying  assumptions  of  previous  controllers 
cannot  be  used,  we  assume  that  bus  work  represents  300  meters  of  cabling.  Cabling 
is  then  paralleled  to  meet  the  MVDC  current  demands  of  the  connected  zone.  We  choose 
300  meters  since  that  distance  represents  the  fore-aft  length  of  an  aircraft  carrier,  the 
U.S.  Navy’s  largest  platform.  We  consider  this  to  be  a  worst-case  assumption,  since 
Figure  10  shows  that  PGMs  and  intermediate  DC-DC  converters  are  distributed 
throughout  the  ship. 

The  input  filter  is  selected  to  be  a  series-RC  filter  connected  from  the  input  of  the 
intennediate  DC-DC  converter  to  ground.  This  filter  was  selected  because  it  adds  a 
stabilizing  pole  to  the  circuit. 

The  intermediate  DC-DC  converter  circuit  models  are  identical  to  Figure  4.  We 
chose  a  fixed  duty  cycle  for  the  intermediate  DC-DC  converters,  as  this  effectively 
renders  them  as  “DC  transformers.”  This  choice  prevents  us  from  simplifying  the 
intennediate  DC-DC  converters  into  ideal  CPL  as  other  authors  have  done.  The  critical 
inductance  and  minimum  capacitance  for  each  intennediate  DC-DC  converter  are 
calculated  assuming  a  1 .0-kHz  switching  speed. 

Finally,  we  consider  the  point-of-load  converters  and  ESD.  The  point-of-load 
converters  are  assumed  to  be  ideal  CPL.  This  assumption  seems  appropriate  since  the 
loads  operating  on  the  zone  buses  operate  at  lower  voltages  than  the  converters  connected 
directly  to  the  MVDC  bus.  This  allows  the  point-of-load  converters  to  operate  at  much 
higher  switching  frequencies  and  exhibit  CPL  behavior.  The  ESD  are  assumed  to  be 
HESS  with  both  fast  dynamics  and  high  power  capacity.  In  order  to  avoid  modeling 
another  complex  system,  the  HESS  is  assumed  to  behave  as  a  controlled  current  source  as 
has  been  done  by  previous  authors.  We  assume  that  HESS  are  connected  to  the  zone 
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buses  by  switching  power  supplies.  Just  as  with  the  CPL,  we  assume  that  the  lower 
voltage  on  the  zone  buses  allows  faster  switching  speeds.  Here  we  have  chosen  16.0-kHz 
switching  speeds  for  the  HESS.  This  number  was  selected  because  it  is  significantly 
faster  than  the  PGM  at  1 .0  kHz  and  because  the  switching  period  divides  into  the  PGM 
switching  period  by  a  convenient  integer  value. 

The  complete  experimental  circuit  model  is  shown  in  Figure  12.  Circuit 
component  names  are  in  black  lettering,  while  state-variables  are  written  in  red  with 
direction  arrows  showing  the  positive  flow  of  currents  and  +/-  symbols  indicating 
polarity  for  voltages.  The  Rgx  and  Lgx  values  indicate  the  PGM  output  impedances.  Since 
all  four  PGMs  are  in  parallel,  the  PGM  output  capacitors  have  been  combined  into  one 
equivalent  bus  capacitor,  Cbus.  The  Rzx,  Lzx,  and  Czx  values  indicate  the  short-wire 
impedance  values.  The  parameters  Rcix  and  C\ix  are  for  the  damper  filter  components.  The 
parameters  Dj,  £>?,  and  D?  are  the  CCM  calculated  duty  cycles  of  the  intennediate  buck 
converters.  The  duty  cycle  values  relate  the  buck  converter  input  current  to  the  buck 
converter  output  voltage.  The  parameters  Lbx  and  Cbx  are  the  inductances  and 
capacitances,  respectively,  for  each  of  the  intermediate  buck  converters.  The  controlled 
current  sources  with  downward  pointing  arrows  are  the  point-of-load  CPL,  while  the 
currents  sources  with  upward  pointing  arrows  are  the  HESS. 

D.  SUMMARY 

In  this  chapter,  we  stepped  through  the  process  of  converting  a  proposed 
NAVSEA  electric  warship  distribution  system  into  an  experimental  circuit  model.  The 
key  features  of  the  NAVSEA  model  were  extracted  and  used  to  make  a  simplified  block 
diagram.  From  the  simplified  block  diagram,  circuit  topologies  and  properties  were 
assigned.  The  final  topology  consists  of  four  average-value  model  DC-DC  converters 
connected  to  a  12.0-kV  MVDC  bus.  Also  connected  to  the  MVDC  bus  are  three  load 
zones.  Each  load  zone  has  cable  impedance,  an  input  RC  filter,  and  an  average  value 
buck  converter  operating  in  CCM  at  fixed  duty  cycle.  On  the  intennediate  buses,  all  loads 
were  assumed  to  be  CPL,  and  the  bus-connected  HESS  was  modeled  as  an  ideal 
controlled  current  source. 
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Figure  12.  Experiment  Circuit  Model 
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IV.  ADAPTIVE,  MULTI-RATE  LINEAR  QUADRATIC 

REGULATOR 


A.  LINEAR  QUADRATIC  REGULATOR 

The  source  document  for  this  section  is  the  Donald  Kirk  text  [36].  LQR  is  an 
optimal  control  technique  for  linear  systems.  One  of  its  attractive  features  is  that  it  can 
be  implemented  recursively  in  an  efficient  manner  using  dynamic  programming.  The 
plant  is  described  in  state-space  form  by  the  system  of  first  order  linear  differential 
equations  as  in 

x(t)  =  A(t)x(t)  +  B(t)u(t) ,  (4.1) 

with  A  representing  the  plant’s  state  matrix,  B  representing  the  feedback  matrix,  x 
representing  the  state -vector,  and  u  representing  the  input  vector.  The  performance 
measure  to  be  minimized,  J,  is  represented  by  the  integral  equation 

1  i  '/ 

J  =  -xT(tf)Hx(tf)  +  -  l  xT  (t)Q(t)x(t)  +  uT  (t)R(t)u(t)~^dt 

2  2  -  ,  (4.2) 

where  H  is  the  boundary  cost  matrix,  Q  represents  the  state-error  cost  matrix,  and  R 
represents  the  input  cost  matrix.  The  matrices  H  and  Q  must  be  real,  symmetric,  and 
positive  semi-definite.  The  matrix  R  must  be  real,  symmetric  and  positive  definite.  This  is 
to  ensure  that  the  inverse  to  R  exists.  For  a  problem  with  finite  time  duration  (tj  <  oo),  the 
optimal  control  is  found  by  recursively  solving  the  Riccati  equation, 

K{t)  =  -Kit) A{t)  -  AT  it)Kit)  -  Qit)  +  Kit)Bit)R  l  it)K(t)  (4  3) 

for  Kit)  starting  at  t/with  Kit f)  =H  and  then  working  backward  in  time  until  to,  with  the 
optimal  input  u(t)  computed  by 

uit)  =  -R\t)BTit)Kit)xit)  (44) 

In  the  special  case  where  certain  conditions  are  met,  Equations  (4.2)-(4.4)  can  be 
dramatically  simplified  while  still  preserving  optimal  state  trajectories.  The  special 
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conditions  are:  1.  the  plant  is  completely  controllable;  2.  A,  B,  Q,  and  R  are  constant 
matrices;  3.  H  =  0;  and  4.  tj  =  oo.  Specifically,  in  the  limit  of  tf  approaching  infinity,  K(t) 
becomes  a  constant  matrix.  For  this  special  case,  called  infinite-horizon  LQR,  Equation 
(4.2)  simplifies  to 

1 00 

J  =  —  j[xT  (t)Qx(t)  +  uT  (t)Ru(t)~^dt 

2  1 0 

while  Equation  (4.3)  becomes 

0  =  -KA  -  AtK  -Q  +  KBR  lK 

and  Equation  (4.4  is  rewritten  as 

u(t)  =  -RlBTKx(t ) 

While  the  completely  controllable  case  provides  optimal  trajectories,  LQR  can 
still  be  used  to  stabilize  the  controllable  modes  in  plants  that  are  not  fully  controllable. 
In  this  vein,  we  expand  the  class  of  problems  for  which  LQR  is  applicable  to 
stabilizable  systems. 

B.  MULTI-RATE  LQR 

The  multi-rate  LQR  controller  theory  presented  in  this  section  is  derived  from 
[37].  Analysis  and  design  of  multi-rate  LQR  controllers  begins  with  the  assumption  that 
the  system  can  be  modeled  as  linear  and  time  invariant  (LTI)  over  a  short  enough  time- 
period  T.  Over  this  short  time-period,  we  may  model  the  plant  with  a  series  of  difference 
equations  as  in 

x[k  +  l\  =  F[k]x[k]  +  G[k]u[k\  for k=  1,2, 3...  ,  (4.8) 

where  x  denotes  the  state  vector,  F  is  a  constant  matrix  that  denotes  the  state  matrix,  G  is 
a  constant  matrix  that  denotes  the  feedback  matrix,  u  is  the  input  vector  and  k  denotes  the 
time  index.  If  the  system  is  periodic  over  a  sequence  of  L  time  steps,  then  a  block-cyclic 
matrix  equation  of  the  form  of 


(4.5) 


(4.6) 


(4.7) 
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Z[i  + 1]  =  AZ[i\  +  BU[i\ 


(4.9) 


where 

Z\i\  —  [Xj ,  xL ,  , ...,  x2  ]  U\i  ]  —  [u,  ,uL,uL_^,...,u2\ 
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characterizes  the  behavior  of  the  system.  For  Equation  (4.9),  the  vector  Z  is  the  super¬ 
state  vector,  A  is  the  super-state  matrix,  B  is  the  super-feedback  matrix,  and  U  is  the 
super-input  vector.  Controllability  and  stabilizability  of  the  Z-step  block-cyclic  system 
may  be  found  by  the  classical  methods. 

Under  this  /.-cyclic  regime,  we  see  that  it  is  possible  to  incorporate  time-varying 
components  such  as  G[k]  in  a  plant  with  inputs  operating  at  multiple  switching 
frequencies.  In  such  a  plant  where  the  switching  events  are  synchronized  and  the 
switching  rates  of  the  various  inputs  may  be  related  by  integer  numbers  of  a  common 
time-step,  the  input  matrix  G[k]  for  this  system  is  periodic  and  this  method  may  be 
employed. 


Computing  block-cyclic  LQR  inputs  is  now  quite  straight-forward.  Following  the 
method  of  Part  A,  we  may  use  the  block-cyclic  A  and  B  matrices  from  Equation  (4.9) 
with  the  block-cyclic  Q  and  R  matrices  given  as 
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to  compute  the  block-cyclic  K  matrix  using  the  discrete-time  Riccati  equation, 

0  =  AtKA  +  Q-AtKB(R  +  BtKBY\AtKB)t  (4 

The  solution  to  the  Riccati  equation  yields  a  block-cyclic  K  whose  components 
Ki-Ki  are  given  as 


Kx  0  0  0  •••  0 

0  K,  0  0  •••  0 

0  0  Kl_  ,  0  •••  0 

:  :  :  :  0 

0  0  0  0  K,  0 

0  0  0  0  0  K2 


(4.12) 


The  block-cyclic  gains  matrix  is  obtained  by  solving  the  block-cyclic  version  of  Equation 
(4.7)  given  by 

G  =  -(R  +  BtKBY'BtKA  13) 

where 

Gx  0  0  0  •••  0" 

0  Gl  0  0  •••  0 

0  0  G,  ,  0  •••  0 

G  = 

:  :  :  :  0 

0  0  0  0  G3  0 

0  0  0  0  0  G2 


The  input  vector  u  calculated  by  the  multiplication  of  the  state -vector  x  by  the  appropriate 
sub-cycle  gain  matrix  Gi  of  the  sequence  as  in 

u[k]  =  Gtx[k]  for  i  =  1,2,  (4. 14) 

Once  the  super  gain  matrix  G  and  all  of  the  G,  matrices  are  found,  the  input  vector 
u  for  time  index  is  found  by  cycling  through  the  sequence  of  G,  matrices  and  multiplying 
them  by  the  state-vector  as  in  Equation  (4.14). 
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This  method,  which  we  refer  to  as  Periodic-Discrete  LQR  (LQR-PD),  has  a  solid 
theoretical  foundation  and  provides  great  value,  especially  when  determining 
controllability.  LQR-PD  is  quite  computationally  costly,  especially  for  high  order 
systems  or  systems  with  a  large  L  number  of  sub-cycles  in  the  sequence.  In  addition  to 
requiring  significant  computation  resources,  we  found  that  the  MATLAB  Riccati  solver 
dare()  had  difficulty  resolving  large-order  systems.  The  dare()  function  issued  warnings 
regarding  “ill-conditioned  matrices”  and  “eigenvalues  too-close  to  the  origin.” 
Furthermore,  for  some  large  systems,  dare()  could  not  resolve  any  solution  at  all. 

Another  drawback  of  LQR-PD  pertains  to  dynamic  systems.  When  the  state 
matrix  A  changes,  such  as  during  a  power  transient  in  an  electrical  distribution  system, 
LQR-PD  is  limited;  it  can  only  make  adjustments  once  every  super-cycle  ( L  indices). 
Depending  on  the  period  of  the  super-cycle,  LQR-PD  may  not  be  able  to  adequately  react 
to  dynamic  events. 

C.  SELECT  MATRIX  LQR 

Select-Matrix  LQR  (LQR-SM)  is  very  similar  to  LQR-PD  except  that  input  gain 
matrices  are  calculated  at  every  sub-cycle  using  sub-cycle  A,  B,  Q,  and  R  matrices. 
Calculating  in  this  fashion  provides  nearly  equivalent  results  to  LQR-PD  while  producing 
advantages  in  computational  complexity  and  microcontroller  memory  requirements. 
LQR-SM  as  an  adaptive,  multi-rate  controller  for  high-order  non-linear  systems  was  first 
presented  in  [38]. 

1.  Description 

To  overcome  the  limitations  of  LQR-PD,  we  need  only  make  two  minor  changes 
to  the  computation  routine.  First,  rather  than  compute  all  of  the  super-cycle  gains  using 
very  large  matrices  as  LQR-PD  does,  we  simply  shift  to  computing  sub-cycle  gains  using 
the  A,  B,  Q,  and  R  matrices  for  the  relevant  sub-cycle.  Second,  rather  than  cycle  through 
several  B  matrices,  we  choose  to  maintain  a  common  B  matrix  and  use  sub-cycle  specific 
R  matrices.  Calculating  gains  on  the  sub-cycle  level  improves  the  speed  of  calculation 
due  to  the  Riccati  solver  (  MATLAB  careQ  or  dare()  )  handling  much  smaller  matrices. 
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This  improvement  in  calculation  speed  is  quite  important  when  feedback  gains  must 
adjust  to  changes  to  the  state  matrix  A  in  real-time,  such  as  dynamic  loading  of  a 
shipboard  electrical  distribution  system. 

R  matrices  are  required  to  be  symmetric  positive  definite;  to  meet  this  requirement 
while  minimizing  the  number  of  design  variables,  we  choose  a  diagonal  R  matrix.  This 
choice  both  simplifies  the  number  of  design  variables  that  must  be  selected  by  the 
designer  but  also  simplifies  the  relationship  between  R  matrix  values  and  regulator 
output.  In  this  way,  the  cost  is  a  function  of  uf,  112  ,  etc.  and  not  of  cross  tenns. 

In  order  to  effect  the  correct  multi-rate  response  from  the  regulator,  each  separate 
sub-cycle  has  an  associated  R  matrix.  Those  inputs  that  are  changing  value  during  the 
switching  sub-cycle  are  assigned  a  diagonal  R  matrix  value  selected  by  the  designer, 
while  those  inputs  that  are  not  changing  during  the  sub-cycle  are  assigned  a  very  large  R 
matrix  value.  For  this  large  R- value,  we  choose  a  value  100-1000  times  larger  than  the 
largest  value  assigned  to  inputs  that  are  changing  during  the  sub-cycle.  With  this  large 
penalty,  the  Riccati  equation  maximizes  use  of  the  unrestricted  inputs  in  order  to  produce 
the  optimal  state  trajectory  while  minimizing  use  of  the  restricted  inputs.  Practice  has 
shown  that  the  restricted  input  values  calculated  in  this  way  are  often  so  small  as  to  be 
nearly  zero.  In  reality,  the  restricted  input  value  is  zero.  Since  the  calculated  value  of  the 
restricted  input  is  a  close  approximation  to  the  actual  value  of  the  input,  there  is  no  great 
deviation  from  reality,  and  the  integrity  of  our  system  is  preserved. 

To  further  generalize  the  application,  we  consider  a  non-linear  system.  The 

analysis  of  Part  B  of  this  section  assumed  that  the  system  could  be  considered  LTI  over  a 

short  time  period  T.  If  the  time  period  T  is  short  enough,  a  non-linear  system  may  be 

approximated  by  a  linear  system.  To  do  so  only  requires  that  the  non-linear  elements  in  A 

and  B  be  approximated  by  versions  of  those  elements  that  are  linearized  about  the 

instantaneous  operating  point.  This  is  where  the  “adaptive”  part  of  the  controller  comes 

in.  Sensors  continuously  measure  the  system  state  variables.  From  the  state  variables,  we 

may  infer  CPL  power  using  the  system  differential  equations,  Ohm’s  law  and  the 

relationship  between  power,  voltage  and  current.  Given  the  measured  CPL  power,  the 

CPL  impedance  is  periodically  linearized.  Using  the  linearized  impedance  of  the  CPL  in 
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place  of  the  actual  non-linear  CPL  impedance  permits  us  to  form  the  linear  differential 
equations  needed  to  perform  LQR.  As  the  CPL  impedances  are  periodically  linearized, 
the  LQR  controller  optimal  feedback  gains  and  input  values  adapt  to  changing  CPL 
impedance. 

2.  Execution  Cycle 

To  better  illustrate  LQR-SM,  we  describe  the  full  execution  cycle  for  one  sub¬ 
cycle.  The  execution  cycles  described  below  are  shown  in  Figure  13,  and  the  controller’s 
block  diagram  is  shown  in  Figure  14. 


1 


Figure  13.  LQR-SM  Execution  Cycle 


Figure  14.  Control  Block  Diagram 


The  first  step  in  the  process  is  to  derive  the  state-space  model  of  the  plant. 
Starting  with  the  physical  model  of  the  plant,  we  derive  the  system  of  first-order 
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differential  equations.  If  the  system  includes  non-linear  differential  equations,  these 
equations  must  be  linearized  about  the  instantaneous  operating  point.  The  system  must 
also  be  fully  observable;  this  requirement  facilitates  both  the  use  of  LQR  as  well  as  the 
calculations  required  to  perfonn  necessary  linearization  of  the  differential  equations. 

Once  the  state-space  model  of  the  system  is  derived,  we  implement  a  change  of 
state-variables.  LQR  regulates  all  state-variables  to  the  zero  condition;  thus,  any  state 
variables  with  a  non-zero  steady-state  value  need  to  be  level-shifted  to  zero. 

Take,  for  example,  the  state  variable  II,  which  is  the  current  through  an  inductor, 
for  an  electrical  system.  The  current  R  should  be  decomposed  into  its  DC  and  small- 
signal  values:  h  =  ii-ss  +  h-DC-  The  state  variable  that  we  want  to  regulate  to  zero  is 
therefore,  we  perform  a  change  of  state  variables  to  define  our  regulated  state  variable  as 
=  //  -  Il-dc.  We  perform  this  change  of  variables  for  all  of  our  state  variables,  making 
appropriate  changes  to  the  A  and  B  matrices.  This  allows  the  designer  to  essentially  deal 
with  the  small-signal  and  DC  portions  of  the  control  problem  separately.  Now,  we  have 
the  A  and  B  matrices  as  well  as  the  state-vector  x  for  the  sub-cycle  of  interest. 

Before  we  can  solve  the  Riccati  equation,  Equation  (4.6),  and  obtain  our  input  via 
Equation  (4.7),  we  must  select  the  appropriate  R  matrix  for  the  switching  cycle.  There 
may  be  up  to  N  /^-matrices  for  a  multi-rate  system  with  N  sub-cycles;  however,  in 
practice  there  are  often  fewer.  Take,  for  example,  a  system  where  one  input  is  updated 
every  sub-cycle  while  a  second  input  value  is  updated  only  every  four  sub-cycles.  With 
four  sub-cycles  per  super-cycle,  we  could  maintain  Rj,  Ri,...,  Rn  matrices,  but  this  is  not 
necessary  since  R2-R4  are  in  fact  identical.  For  a  system  implemented  on  a 
microcontroller,  this  represents  a  memory  savings  over  LQR-PD.  This  advantage  in 
speed  and  memory  extends  even  further  when  we  consider  systems  that  have  large 
numbers  of  sub-cycles  per  super-cycle.  LQR-PD  computational  complexity  grows  as  the 
square  of  the  number  of  sub-cycles.  We  see  that  N  sub-cycles  results  in  an  NxN  matrix  of 
matrices  for  LQR-PD,  while  the  computational  complexity  of  LQR-SM  remains  constant. 
The  only  additional  cost  is  in  the  memory  required  to  store  additional  R  matrices  and  the 
software  to  sequence  those  matrices  appropriately. 
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We  assume  a  constant  Q  matrix  that  is  common  to  all  sub-cycles.  The  selection  of 
appropriate  Q  and  R  matrices  is  left  for  a  later  section. 

Now  that  A,  B ,  Q,  R,  and  x  are  determined  for  the  sub-cycle,  we  input  these  values 
into  Equation  (4.6)  and  obtain  the  small-signal  input  values  u  from  Equation  (4.7).  If  any 
of  our  input  values  also  have  a  DC  component,  the  DC  value  must  be  added  to  the  small- 
signal  value  obtained  from  Equation  (4.7)  in  order  to  obtain  the  total  desired  input  value. 

Now  that  our  input  values  are  known,  we  can  update  the  device  input  values. 
Switching  type  devices,  such  as  DC-DC  converters,  are  typically  governed  by  a  pulse- 
width  modulated  (PWM)  duty  cycle.  For  such  a  case,  we  must  convert  the  desired  device 
output  value  to  the  appropriate  PWM  duty  cycle  for  that  device.  Any  devices  that  are  not 
updated  during  the  computation  cycle  are  assumed  to  hold  a  constant  output  value  or 
constant  duty  cycle. 

3.  Rotating  B-Matrix 

Previously,  we  described  a  method  of  LQR-SM  where  a  common  B  matrix  is  used 
and  the  R -matrix  is  selected  for  the  specific  sub-cycle  combination  of  adjustable  inputs. 
An  alternate  method  of  LQR-SM  is  alike  in  every  way  except  R-matrices  are  developed 
for  each  unique  combination  of  adjustable  inputs  and  the  R -matrix  is  held  common 
among  all  subcycles.  In  this  implementation,  the  R -matrices  developed  are  identical  to 
those  developed  for  LQR-PD  where  the  non-adjustable,  restricted,  inputs  during  the  sub¬ 
cycle  have  their  R-values  set  to  zero. 

The  rotating  R-matrix  implementation  of  LQR-SM  produces  results  much  closer 
to  those  produced  by  LQR-PD  but  suffers  from  reduced  ROAs.  The  exact  mechanism  for 
this  reduction  in  ROA  is  not  presently  known.  We  conjecture  that  rotating  B  LQR-SM 
has  a  more  rigid  structure  owing  to  the  additional  zeroes  in  the  B  matrices.  Unlike  the 
rotating  R  implementation  of  LQR-SM,  the  Riccati  solver  has  less  trade-space  computing 
the  optimal  state  trajectory.  This  increased  rigidity  may  result  in  greater  “brittleness,” 
which  is  manifest  as  smaller  ROAs. 
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D.  COMPARISON  EXAMPLE 

Now  that  we  have  a  basic  description  of  LQR-SM,  we  next  desire  to  compare  its 
performance  to  LQR-PD  and  determine  the  strengths  and  weaknesses  of  each  controller 
relative  to  one  another.  To  do  so,  we  construct  an  example  problem  and  compare  the 
results.  We  then  use  these  results  to  justify  the  conclusion  that  LQR-SM’s  differences  in 
performance  are  outweighed  by  its  advantages  in  computational  complexity.  The 
comparison  performed  in  this  section  was  originally  presented  in  [39]. 

1.  Test  Model 

The  example  we  use  to  illustrate  the  differences  between  LQR-SM  and  LQR-PD 
is  of  a  truncated  version  of  the  NAVSEA  circuit  model  developed  in  Chapter  III.  Figure 
15  is  an  illustration  of  the  truncated  circuit  model.  The  truncated  system  is  composed  of 
two  PGMs.  One  PGM  is  rated  to  deliver  40.0  MW  of  power,  while  the  second  PGM  is 
rated  to  deliver  10.0  MW  of  power.  PGMs  are  modeled  as  controlled  voltage  sources 
with  an  RL  output  impedance.  PGMs  are  connected  directly  to  a  12.0-kV  MVDC  bus. 
The  PGMs  each  operate  at  a  switching  speed  of  1.0  kHz;  thus,  the  PGM  inputs  are 
adjusted  every  1.0  ms. 

The  bus  has  a  buffer  capacitor.  Two  load  zones  are  connected  to  the  MVDC  bus. 
Each  load  zone  is  composed  of  input  cabling  modeled  as  an  RLC  section,  an  RC  input 
fdter,  an  average-value  model  of  a  buck  converter,  a  CPL,  and  a  HESS.  The  CPL  is 
modeled  as  a  controlled  current  source  whose  current  is  equal  to  load  power  divided  by 
the  zone  bus  voltage.  In  this  way,  the  power  delivered  to  the  CPL  is  constant.  The  HESS 
is  modeled  as  a  controlled  current  source.  Both  HESS  are  switched  at  8.0  kHz,  resulting 
in  an  update  period  0.125  ms. 

The  first  load  zone  is  regulated  to  maintain  1.0  kV  with  a  CPL  maximum  loading 
of  20.0  MW.  The  second  load  zone  is  regulated  to  maintain  6.0  kV  with  a  CPL  maximum 
loading  of  28.0  MW.  The  system  is  initially  at  steady-state  with  Zone  #1  loaded  at 
15.0  MW  and  Zone  #2  loaded  at  9.0  MW.  This  represents  the  50%  load  condition.  At 
0.0  seconds,  the  loading  is  instantaneously  stepped  to  20.0  MW  in  Zone  #1  and  28.0  MW 
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in  Zone  #2.  This  represents  the  100%  load  condition.  At  20.0  ms,  the  loading  is 
instantaneously  stepped  back  to  the  50%  load  condition.  The  Q  and  R  matrices  are  tuned 
to  ensure  that  the  voltage  transients  do  not  exceed  ±10%  of  their  steady-state  values. 
Circuit  component  values  are  displayed  in  Table  1. 


Table  1.  Example  Circuit  Component  Values 


RgI 

0.25  0 

4/ 

70.5  pH 

Cdi 

1.7  pF 

Rg2 

0.30  0 

LZ2 

47.0  pH 

Cd2 

2.3  pF 

Lgi 

2.00  mH 

c:, 

2.46  pF 

Lbi 

30.6  pH 

Lg2 

1.80  mH 

C:2 

3.69  pF 

Lb2 

926  pH 

Cbus 

4.0  pF 

Rdl 

10  O 

Cbl 

75  mF 

R:i 

3.30  mO 

Rd2 

10  O 

cb2 

1.25  mF 

R-2 

2.20  mO 

Figure  15.  Two-PGM,  Two-Zone  Experimental  Circuit  Model 

2.  Controller  Implementation 

The  first  step  in  developing  the  controller  is  to  derive  the  differential  equations  for 
the  circuit  in  Figure  15.  These  differential  equations  are  given  as 
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where  Ex  is  the  PGM  voltage  for  machine  “x”;  Igx  is  the  series  inductor  current  for  PGM 
“x”;  Vbus  is  the  MVDC  bus  voltage;  Izx  is  the  line  current  to  zone  “x”;  Vzx  is  the  voltage  at 
the  input  to  buck  converter  for  zone  “x”;  Vdx  is  the  voltage  across  the  damper  capacitor 
for  zone  “x”;  hx  is  the  buck  inductor  current  for  zone  “x”;  V/,x  is  the  voltage  on  the  buck 
filter  capacitor  for  zone  “x”;  Isx  is  the  current  injected  from  HESS  “x”;  Px  is  the  CPL 
power  in  zone  “x,”  and  Dx  is  the  duty  cycle  for  the  buck  converter  for  zone  “x.” 

After  linearization  about  the  instantaneous  operating  point,  the  small-signal  value 
of  CPL  resistance  is  represented  as  the  resistance 
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The  state-space  representation  of  the  small-signal  system  becomes  the  linear  equation 


where 


x(t)  =  Ax(t)  +  Bu(t) . 


(4.17) 
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a.  LQR-PD 

Referring  to  Section  B  of  this  chapter,  to  implement  LQR-PD  we  need  to  take  the 
sub-cycle  matrices  and  form  super-cycle,  block-cyclic  matrices. 

The  first  step  in  this  process  is  to  linearize  the  differential  equations  to  fonn  the 
linear  continuous-time  matrices  A  and  B.  Next,  we  form  the  discrete-time  matrices  F  and 
G  from  the  continuous-time  matrices  A  and  B.  We  do  this  by  recognizing  the  relationship 
between  the  discrete-time  difference  equation  and  the  continuous-time  differential 
equation.  We  convert  the  continuous-time  differential  equation  of  Equation  (4.1)  into  the 
discrete-time  representation  of  Equation  (4.8)  using  the  time-step  between  switching 
events  At  by  the  relationships  of 

x{t)  =  Ax(t)  +  Bu(t) 

F  =  I  +  A-  At 
G  =  B- At 

x[k  + 1]  =  Fx[k]  +  Gu[k ]  (4  18) 

We  now  have  the  F  and  G  matrices  we  need  to  fonn  the  block-cyclic  A  and  B  matrices. 
For  the  system  described  by  Figure  15,  the  voltage  sources,  E}  and  E2,  are  switched  at  a 
rate  of  1.0  kHz  on  alternating  half-cycles.  The  current  sources,  Isi  and  Is2,  are  switched  at 
a  rate  of  8.0  kHz;  thus,  we  can  decompose  the  control  into  eight  0.125-ms  sub-cycles  per 
each  1.0-ms  super-cycle.  During  the  first  sub-cycle,  we  switch  Ej,  Isl  and  Is2.  During  the 
fifth  sub-cycle,  we  switch  E2,  Isl,  and  Is2.  For  sub-cycles  two,  three,  four,  six,  seven,  and 
eight,  we  switch  only  Isl  and  Is2.  Switching  in  this  manner  forms  three  distinct  G 
matrices:  Ga  is  defined  for  the  first  sub-cycle,  G«  is  defined  for  the  fifth  sub-cycle,  and 
Gc  is  defined  for  all  other  sub-cycles.  The  matrix  Ga  has  the  same  structure  as  the  B 
matrix  described  in  Equation  (4.17),  but  Ga2i2  is  zero.  In  a  similar  way,  G«  is  also 
analogous  to  B,  but  Gbu  is  zero,  and  Gc  is  also  like  B,  but  Gci,i  and  Gc2,2  are  zero. 

With  F  and  the  set  of  G  matrices  known,  we  may  form  the  block-cyclic  state  and 
feedback  matrices  A  and  B  according  to  the  process  described  in  Section  B  of  this 
chapter.  The  Q  and  R  matrices  are  formed  with  each  entry  given  by 
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Once  A  ,B  ,Q and R  are  formed.  Equation  (4.11)  can  be  solved  with  MATLAB 
dareQ  and  gain  matrices  calculated  with  Equation  (4.13).  Once  the  gain  matrices  are 
known,  then  the  state  vector  x[k]  can  be  multiplied  by  the  appropriate  gain  matrix  for  the 
sub-cycle  to  find  the  desire  input  values. 

Now  that  the  small-signal  input  values  have  been  found,  we  add  them  to  the 
steady-state  values  for  the  inputs  in  order  to  determine  the  total  commanded  voltage  or 
current  for  the  input  devices.  To  enforce  load  sharing  between  the  two  voltage  sources, 
we  set  their  values  as  given  by 


Ex[k]  =  Vref  +0.80  Pcpli+Pcpli  Rgl+ux[k\ 

Vref 

E2[k\  =  Vref  +  0.20  PcpL\ +  PcpL1  Rgl  +u2[k].  (4.20) 

*  ref 

Isx\k\  =  un[k~\ 

J,2[*]  =  «13[*] 

This  ensures  that  PGM  #1  and  PGM  #2  share  load  according  to  their  ratings  and  that 
MVDC  bus  voltage  is  maintained  at  Vref  for  steady-state  operation.  PGM  #1  provides 
80%  of  total  load  current  at  steady  state,  while  PGM  #2  provides  the  remaining  20%.  The 
HESS  current  sources  have  DC  values  of  zero;  therefore,  the  total  current  demanded  is 
simply  the  small-signal  current  demand. 
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b.  LQR-SM  Rotating  R 

To  implement  LQR-SM  control  by  using  rotating  R-matrices,  we  can  re-use  many 
of  the  elements  found  in  LQR-PD.  The  A  and  B  matrices  of  Equation  (4.17)  remain  the 
same,  as  does  the  Q  matrix  of  Equation  (4.19).  Where  LQR-PD  defined  GA,  Gb  and  Gc, 
we  similarly  define  RA,  Rb  and  Rc;  Ra  is  defined  for  the  first  sub-cycle,  Rb  is  defined  for 
the  fifth  sub-cycle,  and  Rc  is  defined  for  all  other  sub-cycles.  The  RA  matrix  has  the  same 
structure  as  the  R  matrix  described  in  Equation  (4.19),  but  RA2,2  is  set  to  10  .  In  a  similar 
way,  Rb  has  the  form  of  R,  but  Rbu  is  set  to  10  .  Likewise,  Rb  is  has  the  form  of  R,  but 
Rci,i  and  Rcia  are  set  to  10  .  The  full  sequence  of  R-matrices  in  the  super-cycle  is 
displayed  in  Table  2. 

Unlike  LQR-PD,  which  is  only  linearized  every  super-cycle,  LQR-SM  is 
linearized  every  sub-cycle.  After  linearization  about  the  instantaneous  operating  point,  for 
a  given  sub-cycle  we  know  A,  B ,  Q,  and  R.  Equation  (4.6)  is  solved  for  the  sub-cycle  with 
inputs  calculated  by  Equations  (4.7)  and  (4.20). 


Table  2.  LQR-SM  Rotating  R  Super-Cycle 


Super-cycle  Sequence 

Sub-cycle 

1 

2 

3 

4 

5 

6 

7 

8 

R  matrix 

Ra 

Rc 

Rc 

Rc 

Rb 

Rc 

Rc 

Rc 

c.  LQR-SM  Rotating  B 

To  implement  LQR-SM  control  by  using  rotating  B  matrices,  we  define  BA,  Bb, 
and  Be  in  much  the  same  way  as  we  defined  GA,  Gb,  and  Gc  for  LQR-PD.  The  BA  matrix 
has  the  same  structure  as  the  B  matrix  described  in  Equation  (4.17),  but  BA2,2  is  zero.  In  a 
similar  way,  Bb  has  the  form  of  B ,  but  Bbi.i  is  zero.  Likewise,  Be  has  the  form  of  B,  but 
Bcij  and  Bc2,2  are  zero.  The  full  sequence  of  R-matrices  in  the  super-cycle  is  displayed  in 
Table  3.  The  A  matrix  is  the  same  as  in  the  rotating  R  implementation  and  is  linearized 
every  sub-cycle.  The  Q  and  R  matrices  are  identical  to  those  in  Equation  (4.19).  As 
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before,  with  A,  Q,  R  and  the  appropriate  B  for  the  sub-cycle  selected,  we  solve  Equation 
(4.6)  and  calculate  the  inputs  for  the  sub-cycle  by  Equations  (4.7)  and  (4.20). 


Table  3.  LQR-SM  Rotating  B  Super-Cycle 


Super-cycle  Sequence 

Sub-cycle 

1 

2 

3 

4 

5 

6 

7 

8 

B  matrix 

Ba 

Be 

Bc 

Be 

Bb 

Be 

Be 

Be 

3.  Results 

Results  comparing  the  performance  of  LQR-PD  and  both  implementations  of 
LQR-SM  are  described  below. 

a.  LQR-PD  versus  LQR-SM  Rotating  R-Matrix  Implementation 

In  Figure  16a,  PGM  voltages  are  compared  between  the  LQR-PD  and  LQR-SM 
controllers.  The  LQR-PD  controller  employs  a  wider  range  of  PGM  voltages  with  more 
oscillation  than  the  LQR-SM  controller.  Despite  this  wider  dynamic  range  of  PGM 
voltage,  the  currents  through  the  PGM  output  inductors,  shown  in  Figure  16b,  are  very 
similar.  These  greater  voltage  oscillations  manifest  in  greater  variation  in  MVDC  bus 
voltage  and  a  longer  MVDC  voltage  settling  time  for  the  LQR-PD  controller,  as 
illustrated  in  Figure  17. 
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Generator  Voltages 


Figure  16.  PGM  Normalized  Voltages  (a)  and  Currents  (b) 


Figure  17.  MVDC  Bus  Voltage 


When  we  compare  the  voltages  in  each  of  the  load  zones,  shown  in  Figure  18a, 
we  see  that  the  voltages  match  very  closely.  Settling  times  have  no  discernible  difference. 
The  buck  inductor  current  transients  in  both  zones,  shown  in  Figure  18b,  also  match  quite 
closely.  The  LQR-PD  current  transient  has  a  high-frequency  oscillation.  This  oscillation 
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may  be  a  response  by  the  HESS#2  current  to  counter  the  larger  MVDC  bus  oscillations 
produced  by  the  LQR-PD  controller. 


Figure  18.  Normalized  Zone  Buck  Voltages  (a)  and  Currents  (b) 


The  oscillations  exerted  by  the  LQR-PD  controller  are  also  reflected  in  the  HESS 
#2  current  transient  of  Figure  19a.  Despite  what  appears  to  be  greater  control  effort,  the 
LQR-PD  controller  requires  60%  less  HESS  energy  in  Zone  #1  and  8%  less  HESS  energy 
in  Zone  #2  to  stabilize  the  transient.  The  HESS  energy  levels  are  displayed  in  Figure  19b. 

An  analysis  of  regions  of  attraction  yields  more  comparisons.  For  each  of  the 
ROAs  shown  in  Figures  20  and  21,  the  circuit  of  Figure  15  is  at  steady-state,  operating  at 
100%  loading  in  both  zones.  For  Figure  20,  the  Zone  #1  current  and  voltage  were 
disturbed  in  the  region  shown.  The  circuit  and  controller  responses  were  evaluated.  If  the 
circuit  returned  to  steady-state  operation  without  allowing  bus  voltage  to  dip  below 
0.0  V,  without  experiencing  a  Riccati  solution  error,  or  without  exceeding  a  40.0-ms 
settling  time,  the  disturbance  was  evaluated  as  within  the  ROA.  Both  Figures  20  and  21 
show  no  significant  difference  in  ROA  between  the  two  controllers. 
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Seconds 


Figure  19.  HESS  Current  (a)  and  Delivered  Energy  (b) 


Initial  CPL  Amps  xio4 


Figure  20.  Zone  #1  Region  of  Attraction 
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Figure  21.  Zone  #2  Region  of  Attraction 

b.  LQR-PD  versus  LQR-SM  Rotating  B-Matrix  Implementation 

When  we  examine  the  simulation  results  for  LQR-SM  with  Rotating  B  matrices 
versus  LQR-PD,  we  see  that  the  two  controllers  behave  nearly  identically  to  one  another. 
The  PGM  voltages,  presented  in  Figure  22a,  commanded  by  each  controller  are  very 
similar  while  the  PGM  output  inductor  currents,  presented  in  Figure  22b,  are  nearly 
identical. 


a) 


Seconds 


Figure  22.  PGM  Voltages  (a)  and  Currents  (b) 
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This  similarity  of  PGM  voltages  and  currents  extends  further  when  we  examine 
the  MVDC  bus  transients  as  shown  in  Figure  23.  Here  again,  the  differences  between  the 
two  controllers  are  minor.  The  LQR-SM  controller  has  a  slightly  smaller  range  of  MVDC 
bus  voltages  and  a  quicker  settling  time.  The  zone  voltages  and  buck  inductor  currents, 
presented  in  Figure  24,  are  nearly  identical  for  both  controllers.  Inspection  of  the  HESS 
currents  and  delivered  energies  in  Figure  25  shows  that  these  two  controllers  are,  again, 
very  similar  with  LQR-PD;  both  controllers  stabilize  the  transient  with  marginally  more 
HESS  energy. 


Figure  23.  MVDC  Bus  Voltage 


Seconds 


Figure  24.  Normalized  Zone  Buck  Voltages  (a)  and  Currents  (b) 
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Figure  25.  HESS  Currents  (a)  and  Delivered  Energy  (b) 

The  major  difference  between  the  two  controllers  appears  when  we  examine  the 
ROAs  presented  in  Figures  26  and  27.  Both  Figures  26  and  27  we  see  that,  in  the 
investigated  region,  the  ROA  for  LQR-SM  with  Rotating  B  matrices  is  smaller  than  for 
LQR-PD.  Comparing  the  ROAs  in  Figures  20  and  21  to  Figures  25  and  26,  respectively, 
we  notice  that  the  LQR-SM  with  Rotating  B  matrices  has  smaller  ROAs  than  both  LQR- 
PD  and  LQR-SM  with  rotating  R  matrices. 

The  regions  that  are  unstable  for  LQR-SM  with  Rotating  B  but  stable  for  LQR- 
SM  with  Rotating  R  fail  the  40.0-ms  settling  time  test.  As  we  have  seen  in  Figures  23-25, 
the  Rotating  B  control  takes  a  bit  longer  to  settle  than  the  Rotating  R  controller,  which  is 
the  most  likely  reason  why  extreme  disturbances  fail  the  40.0-ms  settling  time  test. 
Allowing  a  longer  settling  time  may  demonstrate  that  these  regions  can  converge  to 
steady-state  values. 
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Figure  26.  Zone  #1  Region  of  Attraction 


Figure  27.  Zone  #2  Region  of  Attraction 


As  mentioned  before,  there  is  a  difference  in  the  computational  complexity  of  the 
competing  controllers.  LQR-PD  requires  solving  the  Riccati  equation  for  very  large 
matrices.  For  a  controller  applied  to  an  NxN  state-space  system  with  L  sub-cycles,  the 
Riccati  solver  must  resolve  a  solution  for  NLxNL  matrices  once  every  L  sub-cycles.  The 
LQR-SM  approach  only  requires  the  Riccati  solver  to  resolve  a  solution  for  NxN  matrices 
but  does  so  every  sub-cycle.  For  the  simulation  trials  presented  in  this  section,  MATLAB 
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is  able  to  complete  the  LQR-PD  simulation  in  about  27  seconds.  For  exactly  the  same 
system,  the  LQR-SM  method  can  produce  simulation  results  in  just  1.9  seconds. 

Attempts  to  expand  this  analysis  to  a  larger  model  of  a  twentieth-order  circuit 
problem,  as  shown  in  Figure  12,  failed.  LQR-PD  could  not  be  implemented  at  the 
twentieth-order  level  due  to  the  MATLAB  Riccati  solver’s  inability  to  resolve  a  solution. 
At  the  thirteenth-order  level,  which  is  presented  above,  the  MATLAB  dare()  function 
produced  warnings  regarding  “ill-conditioned  matrices”  and  “eigenvalues  too  close  to 
the  origin.” 

E.  SUMMARY 

LQR-SM  is  a  multi-input,  multi-rate  LQR-based  control.  LQR-SM  can  be 
implemented  as  either  a  rotating  R  matrix  version  or  a  rotating  B  matrix  version.  Both 
versions  of  LQR-SM  controllers  produce  similar  results  to  LQR-PD  but  do  so  much  more 
quickly.  Due  to  LQR-SM  linearizing  the  state  matrix  more  frequently,  both  LQR-SM 
controllers  provide  slightly  better  regulation  of  the  MVDC  bus  than  LQR-PD  with  faster 
settling  times.  The  Rotating  R  version  of  LQR-SM  has  the  best  MVDC  bus  regulation 
and  fastest  settling  times  of  the  three  controllers  discussed  in  this  section.  LQR-PD  used 
less  HESS  energy  to  stabilize  the  example  transient  than  the  LQR-SM  controllers.  Here, 
we  find  a  trade-off  between  controllers.  LQR-PD  is  computationally  more  complex  but 
produces  a  more  “optimal”  solution  by  using  less  input  energy  to  achieve  the  desired 
regulation.  LQR-SM  is  computationally  less  complex  but  trades  that  simplicity  for  a  less 
optimal  use  of  input  energy. 

Region  of  attraction  analysis  showed  that  LQR-PD  and  LQR-SM  with  rotating  R 
matrices  had  nearly  identical  ROAs  but  that  the  Rotating  B  version  of  LQR-SM  had 
smaller  ROAs. 

The  final  advantage  presented  by  LQR-SM  over  the  classical  multi-rate  LQR  of 
LQR-PD  is  robustness.  Solving  high-order  Riccati  equations  often  proves  difficult, 
especially  when  the  matrices  input  to  the  Riccati  solver  are  as  sparse  as  ours.  By  keeping 
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the  matrices  smaller,  LQR-SM  does  not  experience  the  computational  problems 
encountered  when  using  LQR-PD. 
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V.  STOCHASTIC  SEARCH  DESIGN 


In  this  chapter,  we  investigate  two  different  stochastic  search  algorithms: 
simulated  annealing  and  the  genetic  algorithm.  Stochastic  searches  are  useful  tools  for 
exploring  high-order  search  spaces  and  non-linear  systems.  Here,  stochastic  search  design 
is  used  to  design  an  LQR-SM  controller  of  the  rotating  R-matrix  type  that  meets  certain 
transient  dynamic  specifications  while  minimizing  the  stored  energy  requirement  for  the 
NAVSEA  circuit  model  of  Figure  12. 

LQR  controllers  present  a  challenge  due  to  their  large  number  of  design  variables. 
Each  entry  in  the  Q  and  R  matrices  represents  a  free  variable  that  must  be  assigned  a 
value.  The  requirement  that  Q  is  positive  semi-definite  and  R  is  positive-definite  restricts 
choices  to  some  degree  but  still  leaves  quite  a  bit  of  guesswork  for  the  designer.  A 
stochastic  search  algorithm,  when  properly  constrained,  can  provide  the  necessary 
variable  values.  The  automated  nature  of  a  stochastic  search  allows  the  designer  to  set  the 
algorithm  to  work,  turn  his  or  her  attention  to  other  tasks,  then  return  to  check  the  results 
of  the  search.  This  is  a  much  more  efficient  use  of  human  resources  than  iterative  design. 

In  this  chapter,  we  first  discuss  the  state-space  model  of  the  system  under  study. 
Next,  we  develop  the  performance  criteria  that  are  used  to  evaluate  each  of  our  candidate 
solutions.  Finally,  we  implement  both  simulated  annealing  and  genetic  algorithm 
searches  to  design  the  optimum  LQR-SM  controller  for  the  system. 

A.  SYSTEM  MODEL 

The  system  model  used  is  for  the  average-value  model  developed  in  Chapter  III. 

For  review,  there  are  four  power  generation  modules  (PGM),  two  of  which  are  40.0  MW 

while  the  other  two  are  10.0  MW  for  a  total  of  100.0  MW  of  generating  capacity.  The 

main  MVDC  bus  is  regulated  to  12.0  kV.  There  are  three  load  zones  connected  to  the 

12.0-kV  MVDC  bus.  Each  load  zone  consists  of  a  bus  impedance  derived  from  [27],  a 

series-damped  RC  filter  in  parallel  with  the  medium  voltage  side  of  a  DC-DC  buck 

converter,  the  buck  converter,  an  ideal  CPL,  and  an  ideal  controlled  current-source 

HESS.  The  CPL  and  HESS  are  parallel  connected  to  the  low-voltage  side  of  the  buck 
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converter.  The  three  load  zones  are  a  20.0-MW  capacity  CPL  on  a  1.0-kV  bus,  30.0-MW 
capacity  CPL  on  a  6.0-kV  bus,  and  80.0-MW  capacity  CPL  on  a  10.0-kV  bus.  The  first 
two  zones  have  HESS,  but  the  third  zone  does  not.  Total  load  capacity  exceeds  the  total 
generating  capacity;  however,  overload  conditions  are  not  in  the  scope  of  this  study.  The 
circuit  model  was  developed  in  Chapter  III.B  and  first  illustrated  in  Figure  12;  here,  we 
repeat  the  illustration  as  Figure  28  for  the  reader’s  convenience. 


Figure  28.  Distribution  System  Circuit  Model,  Repeated  from  Figure  12 
The  circuit  model  of  Figure  28  yields  the  differential  equations 


V  s2  J 


V  s4  y 


(5.1a) 
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where  Ex  is  the  PGM  voltage,  Igx  is  the  PGM  inductor  current,  Vi,us  is  the  MVDC  bus 
voltage,  I7X  is  the  line  current  to  zone  “jc,”  Vzx  is  the  voltage  at  the  input  to  buck  converter 
“x,”  Vjx  is  the  damper  capacitor  voltage  for  zone  “x,”  I/)X  is  the  buck  inductor  current  for 
zone  “x”  Vhx  is  the  buck  fdter  capacitor  voltage  for  zone  “x,”  and  Dx  is  the  duty  cycle  for 
the  buck  converter  for  zone  “x.”  The  necessary  state  variables  for  the  LQR  controller  to 
properly  regulate  the  system  are  given  by 
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where  F  is  the  steady-state  total  MVDC  current,  Fre/  is  the  MVDC  bus  reference  voltage 
(12  kV),  Iol  is  the  Zone  #1  steady-state  MVDC  current  (Pcpu/Vrej),  I02  is  the  Zone  #2 
steady  state  MVDC  current  ( PcPL2/Vrej ),  4  is  the  Zone  #3  steady-state  MVDC  current 
( PcPL3/Vref ),  I  obi  is  the  Zone  #1  steady-state  LVDC  current  ( Pcpu/Vrefi ),  Iob2  is  the  Zone 
#2  steady-state  LVDC  current  ( PcPL2/Vref2 ),  Fm  is  the  Zone  #3  steady-state  LVDC  current 
(Pcpu/Vrefi).  The  voltage  Vrep  is  the  Zone  #1  LVDC  reference  voltage  (1.0  kV),  Vrep  is 
the  Zone  #2  LVDC  reference  voltage  (6.0  kV),  and  Vref3  is  the  Zone  #3  LVDC  reference 
voltage  (10.0  kV). 
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B.  EVALUATION  CRITERIA 


Next,  we  develop  the  evaluation  criteria  that  are  the  basis  for  choosing  the 
desirability,  or  fitness,  of  each  candidate  configuration.  Each  candidate  configuration  of 
R,  Q  and  capacitor  values  is  evaluated  based  on  dynamic  response  to  a  presumed  worst- 
case  instantaneous  load  step  from  50%  total  load  power  to  100%  total  load  power.  The 
power  transient  involves  power  in  Zone  #1  stepping  from  15.0  MW  to  20.0  MW,  Zone  #2 
stepping  from  5.0  MW  to  30.0  MW,  and  Zone  #3  stepping  from  28.0  MW  to  46.0  MW. 

The  dynamic  performance  specification  selected  is  that  MVDC  bus  voltage  ( %>U5) 
and  the  three  zone  buck  voltages  (%,/,  V/,2,  and  Vbs)  shall  remain  within  10%  of  steady- 
state  value  during  the  transient.  Given  this  dynamic  performance  specification,  it  is  more 
convenient  to  graph  the  voltage  transients  with  respect  to  the  normalized  voltage  values 
(Vbu/Kef,  V b i/V ref h  etc.)  so  that  the  graph  can  be  quickly  analyzed  to  detennine  if  the 
transient  exceeds  limits.  A  convenient  benefit  is  that  all  of  the  voltage  transients  can  be 
graphed  on  the  same  scale  for  ease  of  comparison  and  presentation. 

Another  specification  placed  on  the  transient  is  that  harmonic  content  not  exceed 
certain  values.  Since  all  PGMs  and  HESS  interface  through  switch-type  power 
electronics,  we  desire  that  the  switching  frequencies  and  certain  harmonics  not  exceed 
-60  dB  down  in  energy  spectral  density  from  the  DC  energy.  This  ensures  relatively  clean 
signal  quality  without  excessive  voltage  ripple.  We  check  for  1.0  kHz  and  4.0  kHz 
harmonics  from  the  PGMs  that  switch  at  1.0  kHz  each  on  staggered  quarter  cycles  as  well 
as  16.0  kHz  from  the  HESS. 

While  the  dynamic  performance  and  harmonic  content  criteria  discussed  above 
are  pass-or-fail  criteria,  our  next  criterion  is  a  continuous  scalar  criterion:  total  system 
stored  energy.  The  stored  energy  in  the  system  is  measured  by  the  energy  stored  in 
capacitors  in  the  steady-state  condition  plus  the  peak  energy  delivered  by  each  HESS 
during  the  design  transient.  The  full  scoring  function  is  described  in  the  MATLAB 
function  fitnessQ  in  Appendix  A.  The  stored  energy  fitness  score  is  given  as 

Fitness  =  ^  7.  CxVc  +  max(\VblIHESSldt)  +  max^\Vb2IHESS2dt'j 

Cbus,Cd\,...,Cb3  ^ 
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(5.3) 


c. 


SEED  CANDIDATE  DEVELOPMENT 


Unlike  the  typical  simulated  annealing  or  genetic  algorithm  processes,  we  do  not 
randomly  choose  the  initial,  or  “seed,”  configuration  of  first-generation  genomes.  Instead, 
we  must  choose  the  seed  configuration  quite  carefully.  Because  our  cost  function  has 
both  discrete  and  continuous  components,  it  is  not  only  possible  but  very  likely  that  a 
randomly  generated  seed  fails  the  discrete  evaluation  criteria.  Subsequent  trials,  which 
are  in  the  neighborhood  of  the  unacceptable  seed,  are  also  very  likely  to  be  unacceptable. 
In  this  way,  the  stochastic  search  algorithm  could  potentially  spend  the  entire  set  of 
iterations  exploring  an  unacceptable  region  of  the  variable-space.  To  prevent  this,  we 
choose  a  seed  that  meets  all  of  our  dynamic  and  hannonic  perfonnance  criteria.  For  our 
particular  case,  we  choose  the  Q  matrix  and  R  matrix  diagonal  terms  to  be  uniformly  “1” 
in  value.  Capacitors  are  made  large  enough  such  that  the  dynamic  performance  criteria 
are  met.  By  using  large  capacitors,  we  may  simultaneously  reduce  overshoot  and 
undershoot  peaks,  filter  unwanted  hannonics,  and  expand  the  sizes  of  ROAs. 

In  order  to  appropriately  size  the  capacitors,  we  employ  the  theory  for  second- 
order  systems  with  CPL  discussed  in  the  Chapter  II.  Detailed  reasoning  and  calculations 
for  developing  the  seed  candidate  are  covered  in  Appendix  B. 

D.  SIMULATED  ANNEALING 

Simulated  annealing  is  a  stochastic  search  algorithm  modeled  after  the 
metallurgical  annealing  process.  In  the  annealing  process,  a  hot  piece  of  metal  is 
repeatedly  heated  and  cooled  to  allow  the  atoms  within  the  alloy  to  fonn  the  crystalline 
structure  with  the  lowest  energy  state.  In  theory,  this  minimum  energy  state  has  inter¬ 
atomic  bonds  with  the  least  amount  of  built-in  stress,  therefore  making  the  crystalline 
structure  as  strong  as  possible.  At  high  temperature,  atoms  are  free  to  bounce  around  and 
explore  several  energy  states.  As  the  system  cools,  the  atoms  become  constrained  to 
lower  and  lower  energy  states  until  temperature  is  so  low  that  all  atoms  are  frozen  in 
place. 
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In  simulated  annealing,  we  define  a  function  to  measure  the  “energy”  of  our 
system  or  process.  Along  with  this  energy  function,  we  also  define  a  global  non-physical 
“temperature”  parameter.  The  simulated  annealing  process  begins  by  randomly  choosing 
the  combination  of  variable  values  to  form  the  initial  energy  state.  Then,  in  successive 
iterations,  the  energy  state  is  pennitted  to  vary  randomly  through  a  neighborhood  about 
the  lowest  observed  energy  state.  As  the  “temperature”  cools,  the  size  of  the 
neighborhood  shrinks.  The  algorithm  is  then  stopped  when  a  certain  energy  level  is 
reached,  a  certain  temperature  is  reached,  or  a  certain  number  of  iterations  have  been 
performed  [40]. 

For  the  implementation  of  simulated  annealing  used  in  this  work,  we  start  with  a 
temperature  of  100.0,  have  a  period  of  constant  temperature  for  1600  iterations,  slowly 
lower  temperature  over  the  next  800  iterations,  then  rapidly  cool  temperature  over  the 
final  800  iterations.  Over  the  total  sum  of  3200  iterations,  the  algorithm  is  given  a  large 
number  of  iterations  of  high  temperature  in  order  to  potentially  escape  a  local  minimum 
in  search  of  the  global  minimum.  Over  the  cooling  period,  the  stochastic  search  is 
concentrated  in  a  smaller  and  smaller  neighborhood  in  order  to  find  the  absolute  lowest 
energy  state  within  the  local  area.  A  diagram  of  the  simulated  annealing  process  is 
featured  in  Figure  29.  In  Figure  29,  we  use  2*  to  denote  the  lowest  energy  candidate,  /  to 
denote  the  neighbor  candidate,  T  for  temperature,  and  J  for  the  cost  of  the  candidate.  In 
this  context,  the  terms  “cost,”  “energy,”  and  “score”  are  interchangeable. 

In  order  to  limit  the  search  space,  we  define  upper  and  lower  bounds  for  each  of 
our  free  variables.  To  generate  the  neighbor  %,  we  choose  to  “mutate”  a  portion  of  the 
design  variables  in  2*  by  a  factor  with  random  normal  distribution  with  a  mean  of  one 
and  variance  equal  to  the  current  temperature  divided  by  the  starting  temperature.  Each 
variable  in  /  is  then  limited  by  the  upper  and  lower  bounds  we  previously  set.  The 
probability  of  any  one  variable  becoming  mutated  is  20%.  This  level  of  mutation  was 
selected  so  that  neighbors  would  be  different  from  the  seed  but  not  radically  different. 
The  MATLAB  code  for  the  neighbor  generating  function,  mutateSAQ,  is  provided  in 
Appendix  A. 
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Figure  29.  Simulated  Annealing  Process  Flow  Chart 


The  cost  performance  results  of  simulated  annealing  optimization  are  displayed  in 
Figure  30.  These  results  are  not  impressive.  Simulated  annealing  yielded  approximately 
two-tenths  of  a  decade  in  stored  energy  reduction  over  3200  iterations  of  the  algorithm. 
This  poor  performance  may  be  due  to  selecting  neighbors  that  are  too  similar  to  the  seed 
value,  or  the  algorithm  may  have  been  trapped  in  a  local  minimum.  Inspection  of  Figure 
30  shows  that  cost  reduction  had  essentially  stalled,  so  further  iterations  would  not  have 
improved  performance.  Changing  the  cooling  profile  of  Figure  31  would  likely  not  have 
much  impact  either,  since  the  local  minimum  appears  to  have  been  found  long  before  the 
cooling  period.  After  several  runs  of  simulated  annealing,  the  performance  of  the  routine 
did  not  produce  results  superior  to  a  few  manual  iterations  by  a  human  operator. 
Dissatisfaction  with  simulated  annealing  led  to  an  exploration  of  genetic  algorithm  as  an 
alternate  optimization  tool. 
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Figure  30.  Simulated  Annealing  Cost  Performance 


Figure  31.  Simulated  Annealing  Temperature  Profile 

E.  GENETIC  ALGORITHM 

Genetic  algorithm  is  a  stochastic  optimization  algorithm  that  mimics  biological 
reproduction  processes.  In  genetic  algorithm,  an  initial  population  is  assembled  from 
randomly  generated  points  in  the  available  design  space.  Each  of  these  points  in  the 
design  space  is  called  a  genome.  Each  genome  is  then  evaluated  according  to  a  fitness 
function,  just  as  we  did  with  simulated  annealing.  In  the  next  step,  some  fraction  of  the 
genomes  in  the  generation  is  culled,  leaving  the  best  performing  genomes  in  the 
generation  to  “breed”  and  produce  the  next  generation.  A  certain  amount  of  mutation  is 
injected  into  the  genomes  of  the  next  generation  in  order  to  produce  new  features  not 
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present  in  the  previous  generation.  Finally,  the  process  is  repeated  until  some  termination 
criterion  is  reached.  As  before,  the  termination  criterion  could  be  a  certain  number  of 
generations  or  a  performance  measure  [41]. 

For  our  implementation  in  [42],  we  chose  a  generation  composed  of  16  genomes. 
We  then  refine  the  population  over  200  generations.  These  numbers  were  chosen  so  that 
genetic  algorithm  would  evaluate  3200  candidate  genomes  just  as  the  simulated 
annealing  algorithm  evaluated  3200  candidates.  Each  genome  is  a  vector  of  Q  and  R 
matrix  diagonals  along  with  the  capacitor  values  we  wish  to  minimize.  This  genome  is 
identical  to  the  one  used  for  simulated  annealing.  The  fitness  function  used  in  the  genetic 
algorithm  optimization  is  Equation  (5.3),  the  same  as  was  used  in  simulated  annealing. 
An  additional  feature  of  genetic  algorithm  is  the  concept  of  The  Elite.  The  Elite  is  the 
genome  with  the  best  observed  fitness.  The  Elite  passes  unchanged  from  one  generation 
to  the  next.  In  this  way,  each  generation  maintains  memory  of  the  best  genome.  The 
genetic  algorithm  block  diagram  is  illustrated  in  Figure  32. 

Just  as  with  simulated  annealing,  we  have  to  be  careful  to  choose  an  initial 
generation  that  has  at  least  one  candidate  pass  all  of  the  discrete  fitness  criteria  for 
dynamics  and  harmonic  content.  One  member  of  the  initial  generation  of  genomes  was 
the  same  exact  seed  used  to  begin  the  simulated  annealing  optimization.  The  remaining 
fifteen  members  of  the  generation  were  genomes  where  approximately  20%  of  individual 
traits  were  selected  by  a  unifonn  random  variable  that  spanned  the  entire  range  for  that 
particular  design  variable.  The  remaining  80%  of  variables  were  set  as  identical  to  the 
seed  genome.  After  each  genome  has  been  scored  by  the  fitness  function,  we  cull  the 
worst  half  of  the  genomes.  The  remaining  best  performing  genomes  are  paired  with  one 
another  and  their  genome  variables  mixed  to  form  children.  Genome  variables  are  chosen 
according  to  a  unifonn  random  variable  for  “mixing.”  The  mixing  vector  M  is  an  Nx\ 
vector  of  ones  and  zeros  where  N  is  the  length  of  the  genome.  The  mixing  vector 
complement  is  M  =  1  -M.  Children  are  created  according  to 


Child,  =  PaM  +  PbM 
Child2  =  PaM  +  PbM  ’ 


(5.4) 
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where  Pa  is  the  first  parent  and  Pb  is  the  second  parent.  In  order  to  create  a  truly  new 
generation,  the  four  genomes  (Pa,  Pb,  and  the  two  children)  are  mutated  according  to  a 
mutating  function.  If  the  elite  genome  is  present  in  the  pairing,  the  elite  genome  is  not 
mutated.  The  mutating  function  used  is  very  similar  to  the  neighbor  generating  function 
used  in  simulated  annealing.  The  major  differences  are  that  the  variance  is  fixed  at  0.5 
rather  than  modulated  by  the  global  temperature  variable. 

Various  combinations  of  population  size  and  number  of  generations  were  tried.  A 
population  of  eight  members  and  400  generations  as  well  as  a  population  of  32  members 
and  100  generations  were  also  tried.  Best  results  were  obtained  with  the  population  of  16 
members  and  200  generations. 

The  genetic  algorithm  cost  perfonnance  chart,  Figure  33,  shows  that  genetic 
algorithm  vastly  outperformed  simulated  annealing.  While  simulated  annealing  only 
achieved  two-tenths  of  a  decade  in  stored  energy  reduction,  genetic  algorithm  achieved 
nearly  two-and-a-half  decades  in  cost  reduction.  Specifically,  the  lowest  system  stored 
energy  requirement  determined  by  simulated  annealing  was  1.176  MJ,  while  the  genetic 
algorithm  was  able  to  find  a  configuration  of  only  343  kJ.  The  genetic  algorithm 
optimized  system  requires  only  one-third  the  stored  energy  of  the  simulated  annealing 
optimized  system. 
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Figure  32.  Genetic  Algorithm  Process  Flow  Chart 


Figure  33.  Genetic  Algorithm  Cost  Performance 
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F.  COMPARISON  TO  STATE-FEEDBACK  LINEARIZATION 

This  section  explores  how  well  the  genetic  algorithm-optimized  LQR-SM 
controller  compares  to  another  controller.  State-feedback  linearization  (LSF)  is  a  non¬ 
linear  control  technique  that  is  relatively  straightforward  to  understand  and  implement. 
Several  examples  of  LSF  controlled  systems  exist  in  the  literature  and  were  explored  in 
Chapter  II  of  this  work.  Reference  [22]  forms  the  foundation  for  the  LSF 
implementations  used  in  this  section. 

To  implement  LSF,  we  first  partition  the  circuit  of  Figure  12  into  second-order 
subsystems.  The  first  subsystem  is  composed  of  the  four  PGMs,  the  MVDC  bus,  and  an 
ideal  CPL.  For  this  first  subsystem,  we  make  two  key  assumptions.  The  first  key 
assumption  is  an  identical  ratio  of  PGM  output  resistance  to  PGM  output  inductance  for 
each  of  the  four  PGMs.  This  assumption  allows  us  to  combine  all  four  PGMs  into  one 
single  virtual  machine.  The  second  key  assumption  is  that  all  three  of  the  load  zones  can 
be  represented  by  a  single  ideal  CPL.  By  doing  so,  we  can  use  the  analysis  of  [22]  to 
design  a  second-order  linear  control  law  for  the  PGMs  to  stabilize  the  MVDC  bus. 

The  remaining  subsystems  are  in  the  load  zones.  By  assuming  that  the  MVDC 
voltage  is  constant,  we  can  reduce  each  of  the  zones  into  an  ideal  voltage  source,  source 
inductance,  filter  capacitor,  and  ideal  CPL.  The  three  zones  are  now  reduced  to  second- 
order  systems.  We  can  then  design  second-order  control  laws  for  the  HESS  currents  to 
stabilize  the  zone  voltages. 

1.  PGM  Subsystem 

For  the  subsystem  consisting  of  the  four  PGMs  and  the  MVDC  bus,  the  load 
zones  are  all  approximated  by  a  single  CPL.  The  equivalent  circuit  model  is  shown  in 
Figure  34. 
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If  we  set  the  condition  that  the  ratio  of  Rgi  to  Lgi  for  all  PGMs  are  equal,  we  can  further 
simplify  Equation  (5.5)  to 
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where  x  is  the  ratio  of  Rgi  to  Lgi  and  Leq  is  the  equivalent  parallel  inductance. 
Suppose  we  set  each  of  the  PGM  voltages  to 


Et  =  Vref  +  S^L.  {Kx  (i Vbus  -  V. )  +  K2Vbus  + 


zP 


C  V  C  V 

bus  bus  bus  bus 


K,.„  . 


2  bus 


(5.7) 


The  parameter  S',  is  a  sharing  or  proportioning  constant  that  allows  each  PGM  to  share 
current  in  proportion  to  its  rating.  The  sum  of  all  sharing  constants  must  add  up  to  one. 
The  K]  and  K2  constant  are  the  gains  of  a  proportional  and  differential  controller.  Notice 
that  Equation  (5.7)  has  both  a  linear  part  and  a  non-linear  part.  The  non-linear  part  of 
Equation  (5.7)  cancels  out  the  non-linear  part  of  Equation  (5.6),  linearizing  the  system. 
By  substituting  Equation  (5.7)  into  Equation  (5.6),  we  get 

G  +(r-K,)Vt„  +(— - - K,)(Vbm-Vr,f)  =  0  (5.8) 

QuMcj 


as  the  transfer  function  for  the  linearized  system.  We  then  choose  the  Kj  and  K2  variables 
to  affect  pole  placement.  Using  the  standard  fonn  of 

5"  +  2£co0s  +  (o\  -  0  (5  9) 

for  a  second-order  linear  differential  equation,  we  set  Kj  and  K2  according  to 
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to  obtain  the  desired  system  dynamics. 

2.  Zone  Subsystems 

To  perform  control  of  the  zones,  we  make  the  simplifying  assumption  that  the 
zone  buck  voltage  is  constant.  This  is  actually  a  fairly  accurate  assumption  considering 
that  we  limit  MVDC  voltage  fluctuation  to  less  than  ±10%.  The  simplified  zone  circuit 
diagram  is  shown  in  Figure  35.  As  before,  we  start  with  the  differential  equations  which 
result  from  the  subsystem  circuit  model,  given  as 

*  (K,~Kefx) 

*bx  j 

Lbx 

*L  =  t!-(4~+4) 

h*  .  (5.11) 


Figure  35.  Zone  Subsystem  Circuit  Model 

Next,  we  take  the  time  derivative  of  the  second  equation  of  Equation  (5.11),  then 
substitute  that  into  the  first  equation  of  Equation  (5.1 1),  to  get 
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If  we  then  posit  that  the  HESS  current  Isx  is  equal  to 

4  =  fa(fa  -4/,)  +  fa  j  (K,  -  v^jut 


(5.13) 
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In  this  case,  Ks  and  K4  are  the  respective  gains  of  a  proportional  and  integral  controller. 
We  then  choose  the  Kj  and  K4  variables  to  affect  pole  placement.  Using  the  standard  form 
of  Equation  (5.9),  we  choose  K$  and  K4  according  to 
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to  obtain  the  desired  dynamic  responses. 

3.  LSF  Parameter  Optimization 

In  order  for  the  LSF  controllers  to  compete  on  even  footing  with  the  LQR-SM 
controller,  the  LSF  design  parameters  were  also  optimized  using  genetic  algorithm.  The 
identical  genetic  algorithm  routine  was  used  to  optimize  the  LSF  controllers  as  well  as 
the  capacitors  we  wish  to  minimize.  The  LSF  controlled  system  has  three  separate 
controllers  that  are  simultaneously  optimized:  one  for  the  PGM  subsystem  and  the  two 
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zones  with  HESS.  Each  controller  needs  to  have  the  damping  factor^  as  well  as  the 
natural  frequency  co0  selected.  That  makes  six  design  variables  for  the  controllers.  The 
seven  capacitors  we  wish  to  minimize  increase  the  total  number  of  design  variables  to 
thirteen.  The  design  variables  are  selected  by  genetic  algorithm  with  a  population  of  16 
and  refined  over  200  generations.  The  fitness  function  Equation  (5.3)  is  used  to  score 
each  configuration. 

4.  Results 

The  results  presented  here  are  for  the  design  transient  stepping  from  50%  power 
to  100%  power  then  back  to  50%  power  again.  The  circuit  model  is  given  by  Figure  12 
with  component  values  given  in  Table  4  and  Table  5.  The  circuit  parameter  values 
obtained  for  the  LSF  based  control  scheme  are  given  by  Table  4.  The  genetic  algorithm 
was  used  to  optimize  for  minimum  stored  energy.  The  same  optimization  was  performed 
for  the  LQR-SM  rotating  R  controller,  for  which  circuit  values  are  listed  in  Table  5.  The 
capacitance  values  differ  between  the  two  controllers,  but  all  other  circuit  parameters 
remain  identical.  The  damping  factors  and  natural  frequencies  used  by  the  LSF 
controllers  are  also  contained  in  Table  4.  Of  particular  note,  the  damping  factors  were 
restricted  to  values  between  0.1  and  1.0,  while  the  natural  frequencies  were  restricted 
from  10.0  Hz  to  8.0  kHz.  From  Table  4,  we  see  that  some  of  these  values  were  optimized 
at  the  limits  of  their  respective  ranges.  Expanding  the  range  of  values  may  provide  better 
results  for  the  LSF  controllers. 


Table  4.  LSF  Component  Values 


Cbus 

333  pF 

czl 

2.46  pF 

Cbl 

75  mF 
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Table  5.  LQR-SM  Component  Values 


Rg] 

0.25  Q. 

R=2 

2.20  mff 

Rcl3 

10.0  Q. 

Rg2 

0.30  Q 

Rz3 

1.30  mO 
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0.24  pF 
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10  Q 
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1.9  mF 

The  voltage  transient  comparison  we  see  in  Figure  36  illustrates  that  the  LSF 
control  scheme  provides  a  very  classical  type  response  with  a  decaying  oscillation.  In 
contrast,  the  LQR-SM  response  has  relatively  few  oscillations.  Except  for  Zone  #2,  both 
controllers  have  very  similar  overshoot  and  undershoot  values  as  well  as  very  similar 
settling  times.  This  is  to  be  expected:  the  genetic  algorithm  selected  for  these  traits.  In 
Figure  36a,  we  see  a  classical  critically-damped  response  from  the  LSF  controller  while 
the  LQR-SM  controller  sharply  peaks  on  undershoot.  The  MVDC  bus  voltage  swells  to 
maximum  overshoot  before  quickly  settling  to  steady-state.  In  Figure  36b,  we  notice  the 
LSF  controlled  system  has  very  little  overshoot.  On  the  other  hand,  LQR-SM  uses  the  full 
overshoot  range.  In  Figure  36c,  both  controllers  provide  an  underdamped  response,  but 
the  LQR-SM  controller  returns  to  steady-state  more  quickly  than  the  LSF  controller.  The 
Zone  #3  voltage  response  in  Figure  36d  is  not  regulated  by  any  controller.  In  both  the 
LQR-SM  and  LSF  cases,  Zone  #3  voltage  mostly  tracks  with  MVDC  bus  voltage 
fluctuations. 


77 


Seconds 


Seconds 


Figure  36.  Optimized  Controller  Voltage  Transients  for  MVDC  Bus  (a), 
Zone  #1  (b),  Zone  #2  (c),  and  Zone  #3  (d) 


The  two  controllers  behave  quite  differently,  especially  in  regard  to  PGM  voltage 
control.  The  LQR-SM  controller,  upon  sensing  the  step-change  in  loading,  responds  by 
rapidly  increasing  PGM  voltages,  shown  in  Figure  37a.  This  anticipatory  surge  in  PGM 
voltage  raises  the  PGM  current  quickly,  allowing  the  MVDC  bus  to  better  maintain 
regulation  while  supplying  the  load  zones.  This  rapid  rise  in  PGM  current  better  supplies 
the  load  current  so  that  less  current  needs  to  be  supplied  by  capacitors.  Less  demand  for 
capacitor  current  allows  capacitors  to  be  sized  smaller  while  still  maintaining  bus 
regulation. 
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In  contrast,  the  LSF  controller  is  more  reactive  than  proactive.  The  LSF  controller 
can  only  respond  to  correct  a  drop  in  MVDC  bus  voltage.  The  LSF  controller,  which  is 
decentralized,  cannot  forestall  voltage  swings  through  anticipatory  shifts  in  PGM  current. 
This  implies  that  the  LSF  controlled  system  requires  larger  capacitors  to  effectively 
regulate  the  MVDC  bus. 

The  LSF  controller  PGM  voltages  are  shown  in  Figure  37b.  Notice  how  the  range 
of  values  of  the  LSF  controlled  PGMs  is  much  narrower  than  the  range  of  values  used  by 
the  LQR-SM  controller  in  Figure  37a.  The  LSF  controller  does  not  require  a  large  gain 
range  for  the  PGM  DC-DC  converters,  while  the  LQR-SM  controller  utilizes  a  much 
wide  gain  range.  This  result  implies  that  control  saturation  is  more  likely  with  the  LQR- 
SM  controller. 
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Figure  37.  PGM  Commanded  Voltage  for  LQR-SM  (a),  and  LSF  (b) 


Now  that  we  have  investigated  the  transient  voltage  response,  let  us  consider 
which  system  is  more  efficient  at  utilizing  stored  energy.  In  Figure  38a,  we  see  the 
transient  curve  for  power  delivered  by  HESS  #1  and  Figure  38b  shows  the  transient  curve 
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for  HESS  #2.  In  both  cases,  the  LSF  controller  requires  greater  peak  energy  from  the 
HESS.  This  appears  to  be  due  to  the  LSF  controller  simply  taking  longer  to  overcome  the 
transient.  The  LQR-SM  controller  more  rapidly  returns  zone  bus  voltages  to  steady-state 
values,  thereby  saving  significant  amounts  of  HESS  energy.  These  points  are  illustrated 
in  Figure  38  and  Figure  39.  The  LQR-SM  controller  HESS  #1  energy  peak  is  about  69% 
of  the  LSF  controller  HESS  #1  energy  peak.  The  LQR-SM  HESS  #2  energy  peak  is  about 
51%  of  that  used  by  the  LSF  controller.  Aside  from  the  energy  savings  in  HESS  #1  and 
HESS  #2  control  effort,  the  capacitors  used  by  the  LQR-SM  controlled  system  are 
smaller  as  well.  Overall,  the  LQR-SM  controlled  system  has  a  stored  energy  requirement 
of  only  61%  of  the  LSF  controlled  system  (339  kJ  for  LQR-SM  versus  555  kJ  for  LSF). 
The  differences  in  performance  may  be  attributable  to  the  fact  that  the  LQR-SM 
controller  is  a  centralized  controller  with  full  system  knowledge  while  the  LSF 
controllers  are  decentralized.  The  LQR-SM  controller  can  synergize  all  of  the  available 
inputs  in  order  to  regulate  the  system.  In  contrast,  each  LSF  controller  is  only  acting  to 
regulate  its  local  bus.  The  lack  of  harmony  between  controllers  contributes  to  their 
inefficiency. 
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Figure  38.  HESS  #1  Transient  Response  for  Current  (a),  Power  (b), 

and  Energy  (c) 
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Figure  39.  HESS  #2  Transient  Response  for  Current  (a),  Power  (b), 

and  Energy  (c) 


G.  SUMMARY 

In  this  chapter,  we  explored  stochastic  search  algorithms  as  a  method  to  design 
LQR-SM  controllers.  Both  simulated  annealing  and  the  genetic  algorithm  were  applied  to 
the  task  of  designing  an  LQR-SM  controller  while  simultaneously  minimizing  total 
system  energy  storage. 

The  simulated  annealing  algorithm  did  not  perform  as  well  as  the  genetic 
algorithm  for  our  example  system.  The  genetic  algorithm,  by  contrast,  provided  a  system 
and  controller  that  met  dynamic  and  harmonic  design  goals  with  much  less  stored  energy 
than  was  required  by  the  simulated  annealing  optimization. 

In  producing  an  optimized  LQR-SM  controller  and  system,  we  wished  to  compare 
the  performance  of  the  LQR-SM  controller  to  a  benchmark  controller.  We  selected  LSF 
as  the  benchmark  due  to  its  popularity  in  the  literature  and  relative  ease  of  understanding 
and  implementation.  In  order  to  produce  the  fairest  possible  comparison,  the  LSF 
controller  was  also  optimized  by  the  genetic  algorithm  using  the  exact  same  population 
size,  number  of  generations,  and  fitness  function.  By  comparing  the  minimum  stored 
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energy  configurations  of  both  LQR-SM  and  LSF,  we  found  that  the  LQR-SM  controlled 
configuration  required  only  61%  of  the  stored  energy  required  by  LSF  to  stabilize  the 
transient  from  50%  power  to  100%  and  back  to  50%  power. 

The  LQR-SM  controller  is  superior  to  the  de-centralized  LSF  controllers.  LQR- 
SM  requires  fewer  simplifications  and  assumptions  than  the  LSF  controller.  LQR-SM,  as 
a  centralized  approach,  is  able  to  harmonize  all  available  inputs  to  provide  optimal 
system-wide  regulation  of  all  buses.  Unlike  LSF,  there  is  no  possibility  of  two  different 
de-centralized  regulators  working  at  cross  purposes  (i.e.,  fighting  controllers).  Overall, 
the  LQR-SM  controller  provides  excellent  bus  regulation  while  requiring  smaller 
capacitors  and  lower  capacity  HESS  as  compared  to  an  LSF  controlled  system.  The 
drawback  of  the  LQR-SM  controller  is  greater  computational  complexity  due  to  the 
requirement  to  solve  the  Riccati  equation  for  large  matrices. 
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VI.  PERFORMANCE  ENHANCEMENT 


One  of  the  weaknesses  of  the  LQR-SM  controller  is  that  solving  the  Riccati 
equation  in  real  time  is  computationally  expensive.  For  non-adaptive  systems,  LQR 
controllers  are  typically  solved  off-line,  and  the  feedback  gain  matrix  is  fixed.  By  shifting 
from  a  fixed  gain  matrix  to  an  adaptive  gain  matrix,  we  introduce  the  computational 
expense  of  solving  the  Riccati  equation  every  switching  cycle.  Using  our  knowledge  of 
the  MVDC  distribution  system  we  have  been  studying,  we  explore  reduced-order 
controllers  and  the  design  trade-offs  imposed  by  them.  Since  our  system  is  not  fully- 
controllable,  we  have  many  uncontrolled  states.  Our  system  is  nevertheless  stabilizable, 
so  these  uncontrolled  states  are  stable.  We  can  develop  controller  models  that  omit  or 
consolidate  the  uncontrolled  states,  then  explore  the  effects  these  reduced-fidelity  models 
have  on  controller  performance.  The  reduced-order  controllers  are  then  applied  to  the  full 
twentieth-order  circuit  model  and  optimized  with  the  genetic  algorithm  as  previously 
described. 

A.  ADAPTIVE  17TH-ORDER  LQR-SM  CONTROLLER 

When  examining  the  results  of  the  LQR-SM  controlled  system,  we  can  see  that 
the  MVDC  bus  voltage,  the  voltage  at  the  input  to  the  zone  buck  converters  (Vzj,  1W,  Vz3), 
and  the  voltages  on  the  filter  capacitors  (Vdi,  Vd2,  Vd3),  are  nearly  identical.  If  we  make 
the  assumption  that  the  filter  damper  resistor  Rjx  is  negligible,  we  can  combine  the 
cabling  parasitic  capacitances  with  the  filter  capacitances.  This  simplification  reduces  the 
order  of  the  controller  from  twentieth  order  to  seventeenth  order.  The  seventeenth-order 
circuit  model  is  shown  in  Figure  40  with  the  regions  of  simplification  circled  in  red. 

Just  as  with  the  adaptive  twentieth-order  LQR-SM  controller,  the  adaptive 
seventeenth-order  LQR-SM  controller  is  optimized  to  meet  our  dynamic  and  hannonic 
specifications  with  minimal  stored  energy  using  identical  methods  to  those  described  for 
adaptive  twentieth-order  LQR-SM. 
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Figure  40.  Seventeenth-Order  Controller  Circuit  Model 


B.  ADAPTIVE  11TH-ORDER  LQR-SM  CONTROLLER 

As  with  the  seventeenth-order  controller,  we  again  note  that  MVDC  voltage  and 
the  zone  buck  input  voltages  are  virtually  identical.  From  this,  we  conclude  that  the 
cabling  parasitic  resistances  and  inductances  ( Rzx  and  Lzx)  are  negligible.  By  ignoring  the 
impedances  of  the  zone  cables,  we  can  further  simplify  the  controller  to  eleventh-order. 
As  with  the  seventeenth-order  controller,  we  also  ignore  the  filter  damper  resistance.  This 
allows  us  to  combine  the  PGM  output  capacitors  with  the  filter  capacitors  to  form  a 
single  equivalent  capacitor.  The  eleventh-order  circuit  model  is  shown  in  Figure  41  with 
the  regions  of  simplification  circled  in  red. 

Just  as  with  the  adaptive  twentieth-order  and  seventeenth-order  LQR-SM 
controllers,  the  adaptive  eleventh-order  LQR-SM  controller  is  optimized  to  meet  our 
dynamic  and  harmonic  specifications  with  minimal  stored  energy  using  identical  methods 
to  those  described  for  adaptive  twentieth-order  LQR-SM. 
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Figure  41.  Eleventh- Order  Controller  Circuit  Model 


C.  NON-AD APTIVE  20TH-ORDER  CONTROLLER 

Having  reduced  the  order  of  the  controllers  to  the  maximum  extent,  we  next 
considered  what  performance  trade-off  might  occur  by  calculating  the  feedback  gains  off¬ 
line.  For  this  experiment,  we  used  a  non-adaptive  LQR-SM  controller.  Rather  than 
estimating  power  in  real  time,  we  linearized  the  system  about  steady-state  100%  power 
operation  and  calculated  the  feedback  gain  matrices  for  each  switching  combination.  In 
our  case,  we  have  four  PGMs  and  two  HESS.  Just  as  with  the  adaptive  LQR-SM 
controllers,  we  have  16  sub-cycles  in  a  super-cycle  and  five  different  switching  input 
combinations.  In  the  first  sub-cycle,  PGM  #1  is  switching  along  with  HESS  #1  and  HESS 
#2.  The  next  three  sub-cycles  have  just  the  two  HESS  updating.  The  fifth  sub-cycle  has 
PGM  #2  and  the  two  HESS  switching  followed  by  three  sub-cycles  of  just  the  HESS 
switching,  and  so  on.  This  results  in  five  distinct  gain  matrices.  One  gain  matrix  is  used 
when  switching  PGM  #1  and  the  two  HESS,  the  second  for  switching  PGM  #2  and  the 
two  HESS,  the  third  for  switching  PGM  #3  and  the  two  HESS,  the  fourth  for  switching 
PGM  #4  and  both  HESS,  and  fifth  gain  matrix  is  for  switching  the  two  HESS  and  no 
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PGMS.  Once  the  gain  matrices  are  computed  offline,  they  are  stored  in  memory  and 
called  upon  during  the  appropriate  sub-cycle  to  compute  the  PGM  and  HESS  input 
values. 

Just  as  with  the  adaptive  LQR-SM  controller,  the  non-adaptive  LQR-SM 
controller  was  optimized  to  meet  our  dynamic  and  hannonic  specifications  with  minimal 
stored  energy  using  identical  methods  to  those  described  for  adaptive  twentieth-order 
LQR-SM. 

D.  RESULTS 

1.  Minimum  Energy  Controllers 

a.  Calculation  Speed 

To  compare  the  computational  load  of  each  controller,  the  different  control 
strategies  were  implemented  in  MATLAB  simulations.  While  inputs  to  the  PGMs  and 
HESS  were  recalculated  every  62.5  microseconds,  the  simulation  updated  the  differential 
equations  every  0.1  microseconds.  The  figures  displayed  are  produced  by  sampling  the 
inputs  and  state  variables  at  40.0  kHz.  The  data  streams  were  down-sampled  due  to 
limitations  in  computer  memory.  The  time  duration  of  the  simulations  was  recorded  using 
the  MATLAB  tic  and  toe  functions. 

The  different  controllers’  speed  perfonnance  improved  by  reducing  the  order  of 
the  controller,  as  expected.  The  non-adaptive  controller  was  obviously  the  fastest, 
completing  the  0.1-s  simulation  trial  in  12.8-s.  Since  the  non-adaptive  controller  only 
computes  the  Riccati  equation  once  off-line,  there  are  no  Riccati  equation  computations 
included  in  that  12.8-s  period.  We  can  regard  this  simulation  time  as  the  minimum 
amount  of  time  that  our  MATLAB  script  can  run  on  the  available  desktop  computers. 
The  next  fastest  controller  was  the  eleventh-order  controller,  which  clocked  a  13.6-s 
simulation  run.  This  0.8-s  difference  represents  1600  Riccati  equation  computations  for 
an  average  of  500  microseconds  of  computing  time  for  each  eleventh-order  Riccati 
solution.  The  seventeenth-order  controller  was  the  next  fastest  at  13.7-s.  The  jump  from 
eleventh-order  to  seventeenth-order  only  added  an  additional  0.1  s  of  computing  time. 
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This  averages  to  563  microseconds  per  Riccati  solution.  The  twentieth-order  controller 
was,  as  expected,  the  slowest  at  14.2  s.  This  averages  to  875  microseconds  per  Riccati 
solution.  Calculation  results  are  displayed  in  Table  6. 


Table  6.  LQR-SM  Controller  Computation  Times 


Controller 

Non-adaptive 

1  lth-Order 

17th'Order 

20th-Order 

Seconds 

12.8 

13.6 

13.7 

14.2 

b.  Energy  Storage  Requirement 

The  next  basis  for  comparison  is  determining  if  there  is  a  performance  trade-off 
with  respect  to  energy-storage  requirements  with  reduced-order  controllers.  The 
hypothesis  going  into  this  comparison  was  that  reducing  the  fidelity  of  the  controller 
models  would  result  in  reduced  “optimality”  from  the  controllers;  therefore,  each 
reduction  in  fidelity  of  the  controllers  would  result  in  an  energy  storage  cost  penalty.  In 
Figure  42,  we  compare  the  genetic  algorithm  cost  performance  of  each  of  the  different 
types  of  controllers  through  200  generations  of  optimization. 

The  LSF  controller  is  the  worst  performer  of  the  group  and  is  included  for 
the  sake  of  having  a  baseline.  The  LSF  controller  had  an  energy  storage  requirement  of 
554.9  kJ.  The  three  adaptive  controllers  had  almost  no  differences  in  their  performances. 
The  eleventh-order,  seventeenth-order,  and  twentieth-order  controllers  produced 
minimum  energy  storage  requirements  of  331.5  kJ,  337.0  kJ,  and  331.6  kJ,  respectively. 
The  hypothesis  that  a  loss  of  controller  fidelity  would  result  in  reduction  in  optimality  is 
false.  The  reason  for  the  parity  in  performance  between  the  adaptive  controllers  is  most 
attributable  to  the  fact  that  the  states  that  were  eliminated  from  the  controller  model  were 
uncontrollable  states.  Although  the  voltages  and  currents  eliminated  from  the  controller 
model  may  have  a  minor  effect  on  the  remaining  states,  an  examination  of  graphs  of  V]nts, 
Vzi,  VZ2,  V:j,  Vdi,  Vd2,  and  Vds  shows  that  these  voltages  are  nearly  identical.  Eliminating 
impedances  that  had  very  minor  effects  on  the  circuit  resulted  in  only  minor  effects  on 
the  controller. 
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The  big  surprise  from  this  comparison  is  that  the  non-adaptive  controller  was 
optimized  to  a  minimum  energy  storage  requirement  of  only  306.9  kJ.  The  non-adaptive 
controller  may  have  the  stored  energy  advantage  for  a  few  reasons.  The  first  reason  is  that 
the  non-adaptive  controller  is  designed  for  1 00%  power,  so  it  does  not  have  to  delay  by  a 
super-cycle  before  all  of  the  inputs  are  adjusted  for  the  power  jump.  Step  changes  in 
loading  below  100%  power  level  may  yield  a  different  outcome.  The  second  reason  could 
be  that  the  genetic  algorithm  located  a  more  optimal  result  for  the  non-adaptive 
controller.  The  genetic  algorithm  has  a  great  deal  of  randomness  built  into  it.  This 
particular  optimization  of  the  non-adaptive  controller  may  have  just  been  lucky  with 
regard  to  receiving  the  right  randomizations  and  mutations  to  produce  an  excellent  result. 
Our  results  with  the  genetic  algorithm  have  shown  that  all  trials  produce  roughly 
equivalent  scores,  but  there  is  no  convergence  to  a  single  set  of  parameters.  The  genetic 
algorithm  gives  us  many  “very  good”  solutions  but  does  not  guarantee  the  “best” 
solution. 


Figure  42.  Comparison  of  Controller  Genetic  Algorithm  Fitness 
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c. 


Region  of  Attraction 


The  ROAs  were  calculated  using  the  same  search  technique  used  in  Chapter  IV. 
The  circuits  were  established  at  100%  power  and  the  zonal  voltages  were  disturbed  in  the 
zonal  I-V  planes.  If  the  circuit  returned  to  steady-state  operation  without  allowing  bus 
voltage  to  dip  below  0.0  V,  without  experiencing  a  Riccati  solution  error,  and  without 
exceeding  a  40.0-ms  settling  time,  the  disturbance  was  evaluated  as  within  the  ROA.  The 
goal  is  to  understand  what  trade-offs  can  be  expected  between  the  different  controllers. 

In  the  minimum  energy  configuration,  each  circuit  is  optimized  for  minimum 
stored  energy  requirement.  This  means  that  each  of  these  different  controllers  has  unique 
capacitance  values  for  C}nis,  Cdi,  Cdi,  Cd3,  Cm,  Cb2,  and  C«.  The  differences  in 
capacitance  values  may  play  a  factor  in  the  size  of  the  ROAs.  The  different  capacitance 
values  for  each  minimum  energy  configuration  are  recorded  in  Table  7. 


Table  7.  Minimum  Energy  Configuration  Capacitance  Values 


Capacitor 

Cbus 

cdI 

Cd2 

cdi 

czl 

cz2 

Q, 

cbI 

CM 

CM 

20th 

5.4  |o.F 

0.1  (tF 

0.3  |uF 

20  |iF 

2.5  (uF 

3.7  jaF 

6.2  (tF 

78  mF 

0.7  mF 

2.2  mF 

17th 

3.5  |uF 

18  |iF 

7.8  |iF 

8.0  jaF 

2.5  (uF 

3.7  jaF 

6.2  |iF 

75  mF 

0.7  mF 

2.2  mF 

11th 

3.5  (uF 

4.6  (uF 

2.6  jaF 

0.2  |iF 

2.5  (uF 

3.7  jaF 

6.2  (tF 

75  mF 

0.8  mF 

2.2  mF 

NA 

3.2  |iF 

0.1  |iF 

45  |iF 

1.1  |iF 

2.5  (uF 

3.7  (aF 

6.2  |iF 

75  mF 

0.7  mF 

1.4  mF 

LSF 

333  (uF 

20  |iF 

0.2  |iF 

13  |iF 

2.5  (uF 

3.7  jaF 

6.2  (tF 

75  mF 

0.7  mF 

2.6  mF 

The  Zone  #1  ROAs  are  displayed  in  Figure  43.  Immediately,  we  see  that  the 
ROAs  are  quite  different  from  one  another.  The  three  adaptive  controllers  have  somewhat 
similar  ROAs,  with  the  twentieth-order  adaptive  controller  having  the  smallest  ROA. 
Again,  the  non-adaptive  twentieth-order  LQR-SM  controller  surprises  us  by  having  the 
largest  ROA.  Examining  the  Zone  #1  buck  filter  capacitors  for  each  configuration,  we 
can  see  that  all  of  the  configurations  have  very  similar  Cm  values.  All  of  the 
configurations  use  a  Cbi  of  75  mF  except  for  the  twentieth-order  adaptive  LQR-SM 
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controller,  which  uses  78  mF.  Despite  each  of  these  configurations  having  nearly 
identical  Cbi  values,  the  Zone  #1  ROAs  are  quite  different.  The  differences  may  lie  in 
factors  such  as  controller  feedback  gains  or  systemic  interactions. 

Examining  the  Zone  #2  ROA,  we  have  a  similar  story  as  with  Zone  #1.  Prior  to 
examining  each  of  the  ROAs,  one  would  have  suspected  that  the  higher  fidelity 
controllers  would  have  an  advantage  over  the  lower  fidelity  controllers  or  that  larger  zone 
buck  capacitors  had  an  advantage  over  smaller  capacitors,  but  neither  of  these  predictions 
are  supported  by  the  evidence.  In  the  Zone  #2  ROAs  displayed  in  Figure  44,  the  three 
adaptive  controllers  each  perform  very  similarly,  while  the  non-adaptive  controller  has  a 
noticeably  larger  ROA.  Again,  looking  at  the  Ct2  values  from  Table  7,  we  see  that 
capacitor  size  had  no  discernable  effect  on  ROA  size.  The  controller  action  must  be 
playing  a  significant  role.  Even  though  all  controllers  have  been  selected  through  the 
genetic  algorithm  for  the  same  dynamics,  the  differences  in  feedback  gains  seem  to  have 
a  profound  effect  on  ROAs. 


Figure  43.  Zone  #1  Minimum  Energy  ROAs 


90 


12000 


10000 


w 

o 

> 


8000 


6000 


co 

“  4000 


17th 

11th 

NA 


2000 


0  - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 

0  1000  2000  3000  4000  5000  6000  7000  8000  9000  10000 

Initial  CPL  Amps 


Figure  44.  Zone  #2  Minimum  Energy  ROAs 


The  Zone  #3  ROAs,  shown  in  Figure  45,  are  all  quite  similar.  We  expect  that  Cbs 
size  will  have  more  effect  on  Zone  #3  than  Cbi  or  Cb2  had  on  their  respective  zones  since 
Zone  #3  does  not  have  a  HESS  current  source  to  influence  dynamics.  Despite  all  three  of 
the  adaptive  controllers  having  identical  C«  values,  there  is  a  noticeable,  but  small, 
spread  on  Zone  #3  ROAs.  Even  the  non-adaptive  controller,  with  its  significantly  smaller 
Cb3  value,  still  has  an  ROA  on  par  with  the  three  adaptive  controllers. 


Initial  CPL  Amps 

Figure  45.  Zone  #3  Minimum  Energy  ROAs 
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The  LQR  controller  Q  and  R  matrix  diagonal  values  are  recorded  in  Tables  8  and 
9,  respectively,  for  completeness.  There  is  actually  very  little  difference  between  the 
values  in  Table  9  for  each  controller.  The  differences  in  Q  matrices  are  more  varied,  but 
in  general,  values  are  similar.  It  is  interesting  to  note  that  each  of  these  converters  has  a 
different  solution  for  achieving  the  same  result.  All  four  of  the  LQR-SM  controllers  have 
been  optimized  by  the  genetic  algorithm  for  200  generations  and  each  achieves  a  very 
similar  final  score.  The  genetic  algorithm  performance  curves  of  Figure  42  might  lead 
one  to  believe  that  these  controllers  have  reached  a  global  minimum  or  near  global 
minimum.  We  have  seen  in  actuality  that  each  of  the  genetic  algorithm  routines  has 
converged  not  on  the  unique  “best  solution”  but  on  a  set  of  “very  good”  solutions. 
Despite  all  of  the  solutions  in  the  set  of  “very  good”  solutions  meeting  our  minimum 
energy  performance  goals,  there  can  be  marked  differences  in  the  size  and  shape  of 
ROAs  created  by  each  solution.  This  is  another  potential  area  for  optimization. 
Unfortunately,  generating  a  ROA  takes  approximately  four  hours  of  computation  time, 
making  optimization  in  this  domain  a  costly  endeavor. 


Table  8.  Minimum  Energy  LQR  Controller  Q-values 


Q-Penalty 

20th 

17th 

11th 

NA 

Q-Penalty 

20th 

17th 

11th 

NA 

hi 

1 

l.l 

1 

15 

V,3 

1.2 

1.5 

- 

1.6 

4 2 

1 

1.4 

1 

1.5 

vdI 

1.1 

- 

- 

2.9 

hi 

1 

1.1 

1 

15 

vd2 

2.2 

- 

- 

2.2 

4  4 

1 

1.4. 

1 

1.5 

Vd3 

1.4 

- 

- 

1.1 

V bus 

1 

1.6 

2.2 

1.2 

hi 

1 

3.1 

2.4 

2.9 

hi 

l.l 

6.8 

- 

1.3 

Ib2 

1.7 

4.1 

3.3 

7.4 

1,2 

4.9 

6.7 

- 

1.7 

Ib3 

5.9 

6.1 

1 

1.6 

Iz3 

1.9 

1.9 

- 

24 

vbJ 

2.4 

12 

62 

1 

vzl 

1 

3.3 

- 

2.9 

vb2 

137 

258 

151 

382 

V,2 

1 

2.6 

- 

1.2 

vb3 

1 

14 

1.2 

2.3 
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Table  9.  Minimum  Energy  LQR  Controller  R-values 


Input  R-Penalty 

20th 

17th 

11th 

NA 

PGM  #1/3 

4 

3 

2.4 

6.8 

PGM  #2/4 

1 

2 

1 

1 

Unavailable 

1000 

1000 

1000 

1000 

HESS#1 

0.1 

0.2 

0.7 

0.1 

HESS#2 

1 

1 

1 

1 

2.  Identical  Configuration  Controllers 

Since  the  results  of  the  minimum  energy  optimized  controllers’  ROAs  variance 
seemed  highly  dependent  on  controller  design  parameters,  in  this  section  we  investigate 
differences  in  controller  performance  when  all  controllers  are  using  the  same  capacitor 
values  and  controller  values.  For  this  section,  all  four  LQR-SM  controller  variants  were 
optimized  together  using  200  generations  of  the  genetic  algorithm.  In  this  optimization, 
all  four  controller  variants  were  required  to  use  the  same  capacitor  values.  To  make  the 
comparison  even  closer,  each  controller  uses  identical  Q  and  R  penalty  values  for  state 
variables  and  inputs.  For  example,  for  all  four  controllers,  the  ^-matrix  penalty  on  Vbi  is 
identical.  Even  though  the  reduced-order  controllers  do  not  account  for  some  state 
variables,  at  least  there  is  parity  between  the  controllers  for  all  of  the  state  variables  they 
have  in  common. 

Since  nothing  has  changed  with  respect  to  computation  of  the  respective 
controllers,  the  computation  speeds  of  the  identical  configuration  controllers  are  the  same 
as  those  presented  in  the  minimum  energy  configuration  section. 

The  energy  requirement  for  the  identically  configured  controllers  is  in  line  with 

the  results  of  the  minimum  energy  configurations.  The  three  adaptive  controllers  all 

require  about  the  same  amount  of  total  stored  energy,  while  the  non-adaptive  controller 

requires  somewhat  less  energy.  The  twentieth-order  controller  required  345  kJ,  the 
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seventeenth-order  controller  required  340  kJ,  the  eleventh-order  controller  needed  341  kJ, 
and  the  non-adaptive  controller  required  only  329  kJ. 

When  examining  the  ROAs  of  the  identically  configured  controllers,  the  three 
adaptive  controllers  had  ROAs  that  almost  perfectly  overlapped.  This  fact  is  a  major 
contrast  to  the  amount  of  variability  we  saw  in  Figure  43.  The  non-adaptive  controller  has 
a  more  limiting  ROA,  as  seen  in  Figure  46. 


Figure  46.  Zone  #1  Identical  Configuration  ROAs 

In  Zone  #2,  we  again  see  that  the  adaptive  controllers  all  match  quite  closely  with 
one  another.  The  non-adaptive  controller  actually  has  a  slightly  larger  ROA  in  this  zone, 
though  the  difference  is  much  less  pronounced  than  in  Zone  #1.  The  Zone  #2  ROAs  for 
identical  configuration  are  displayed  in  Figure  47.  In  Zone  #3,  the  ROAs  are  all  quite 
similar,  though  the  non-adaptive  controller’s  ROA  is  noticeably  smaller  than  the  adaptive 
controllers.  The  identical  configuration  Zone  #3  ROAs  are  shown  in  Figure  48. 

From  these  results,  we  can  see  that  with  identical  Q  and  R  matrix  penalty  values, 
the  reduced-order  controllers  behave  almost  identically  to  one  another  except  for 
improvements  in  computation  speed.  There  is  no  measureable  advantage  to  the  higher- 
order  adaptive  controllers.  The  non-adaptive  controller  has  shown,  at  least  in  the 

94 


realizations  explored  here,  advantages  in  computation  speed  and  stored  energy 
requirement.  The  non-adaptive  ROA  size  is  an  area  of  interest.  In  some  cases,  the  non- 
adaptive  controller  has  shown  a  larger  ROA  than  the  adaptive  controllers  but  in  other 
cases  has  shown  a  significantly  smaller  ROA. 


Initial  Amps 

Figure  47.  Zone  #2  Identical  Configuration  ROAs 


Initial  Amps 


Figure  48.  Zone  #3  Identical  Configuration  ROAs 
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Identical  configuration  circuit  component  values  are  enumerated  in  Table  10.  The 
^-matrix  values  are  listed  in  Table  1 1  with  R  values  in  Table  12. 


Table  10.  Identical  Configuration  LQR  Component  Values 


Rgi 

0.25  n 

Rz2 

2.20  mfi 

Rd3 

ion 

Rg2 

0.30  Q 

R-.3 

1.30  mfi 

cdl 

2.73  pF 

Rg3 

0.26  Q 

Lzi 

70.5  pH 

Cd2 

38.2  pF 

Rg4 

0.32  n 

Lz2 

47.0  pH 

Cd3 

17.4  pF 

Lgi 

2.00  mH 

Lz3 

28.2  pH 

Lbi 

30.6  pH 

Lg2 

1.80  mH 

Czl 

2.46  pF 

Lb2 

1.8  mH 

Lg3 

1.95  mH 

cz2 

3.69  pF 

Lb3 

298  pH 

Lg4 

1.71  mH 

Cz3 

6.15  pF 

Cm 

75  mF 

Cbus 

17.2  pF 

Rdl 

10  Q 

Cb2 

0.7  mF 

Rzi 

3.30  m£2 

Rd2 

10  Q 

Cb3 

2.0  mF 

Table  1 1.  Identical  Configuration  LQR  Controller  Q- Values 


Qigi 

l.i 

Qhi 

1.3 

Qvz3 

1.4 

Qlb2 

5.9 

Qlg2 

2.1 

Qlz2 

2.8 

Qvdl 

6.3 

Qlb3 

6.1 

Q& 

1.1 

Qlz3 

2.9 

Qvd2 

8.6 

Qvbi 

4.8 

Qlg4 

2.1 

Qvzl 

2.1 

Qvd3 

2.4 

Q.Vb2 

364 

Qvbus 

14.5 

Qvz2 

1.4 

Qlbl 

5.9 

Qvb3 

1.3 

Table  12.  Identical  Configuration  LQR  Controller  R-Values 


PGM  #1/3 

PGM  #2/4 

ESD  #1 

ESD  #2 

Rmax 

4.4 

1.9 

0.9 

1.0 

1000 

E.  SUMMARY 

In  this  chapter,  we  explored  various  design  trade-offs.  The  MVDC  distribution 
system  we  model  has  several  state  variables  that  very  closely  track  one  another.  These 
state  variables  arise  due  to  modeling  of  bus  impedances  and  input  filters.  If  we  assume 
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that  input  filter  resistance  values  are  negligible,  we  achieve  a  reduced-order  model  with 
only  17  state  variables.  If  we  also  assume  that  bus  impedances  may  be  ignored,  we  can 
eliminate  six  more  state  variables  to  achieve  an  eleven  state-variable  model. 

By  reducing  the  order  of  our  model  through  elimination  of  uncontrollable  modes, 
we  can  reduce  the  size  of  the  A,  B,  Q  and  R  matrices  that  are  input  to  the  Riccati  equation. 
The  reduction  in  matrix  size  significantly  lightens  the  computational  load  of  the  Riccati 
solver.  Computation  time  can  be  further  minimized  by  choosing  a  non-adaptive  LQR-SM 
controller  which  solves  the  Riccati  equation  off-line.  The  excellent  performance  of  the 
non-adaptive  controller  is  aided  by  our  tight  voltage  restrictions.  If  voltage  undershoot 
and  overshoot  were  allowed  greater  range,  the  results  might  have  been  different. 

After  optimizing  each  of  the  controllers  through  identical  genetic  algorithm 
routines,  we  found  no  degradation  in  the  minimum  stored  energy  required  by  each 
controller.  All  four  of  the  controllers  (twentieth-order  model,  seventeenth-order  model, 
eleventh-order  model,  and  non-adaptive  twentieth-order  model)  achieved  roughly  the 
same  stored  energy  value. 

Examination  of  the  ROAs  illustrated  that  ROA  size  is  highly  variable.  Although 
the  minimum  energy  configuration  controllers  had  similar  capacitor  values  for  Q,/,  C/,2, 
and  Ch3,  their  ROAs  in  Zone  #1,  Zone  #2,  and  Zone  #3  were  dissimilar  to  one  another, 
respectively.  These  differences  in  ROAs  illustrate  that  although  all  controllers  perform 
similarly  according  to  the  minimum  energy  fitness  function  Equation  (5.3),  they  do  not 
perform  similarly  with  regard  to  ROA  size.  When  the  controllers  were  constrained  to 
have  identical  penalties  for  like  state  variables,  the  adaptive  controllers  produced  nearly 
identical  ROAs;  however,  the  non-adaptive  controller  ROAs  were  different  from  the 
adaptive  controller  ROAs.  Sometimes  the  non-adaptive  controller  ROAs  were  larger, 
other  times  smaller.  If  there  are  specific  ROA  requirements,  those  requirements  need  to 
be  built  into  the  genetic  algorithm  fitness  function. 
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VII.  CONCLUSIONS 


In  this  dissertation,  we  set  out  to  find  a  suitable  controller  to  regulate  the  voltage 
of  a  hypothetical  naval  MVDC  shipboard  electrical  distribution  system.  That  distribution 
system  would  include  a  zonal  architecture,  CPL,  and  HESS. 

Background  work  demonstrated  how  constant  power  loads  destabilize  electrical 
systems.  CPL  introduce  negative,  non-linear  impedance.  Linearization  of  CPL  impedance 
has  shown  that  care  must  be  taken  in  the  selection  of  circuit  components  in  order  to  avoid 
unstable  eigenvalues.  Non-linear  analysis  demonstrated  that  even  systems  that  are  small- 
signal  stable  may  not  be  globally  stable.  Rather,  they  exist  within  ROAs.  ROAs  can  be 
very  sensitive  to  circuit  component  selection  and  CPL  loading. 

A  review  of  related  work  found  that  many  authors  have  addressed  the  problem  of 
electrical  distribution  systems  with  constant  power  loads  but  with  a  very  limited  focus. 
The  CPL  stabilization  methods  and  controllers  in  the  literature  all  address  the  problem  as 
a  single-input  control  problem,  often  simplifying  the  distribution  system  into  a  simple 
second-order  equivalent  circuit. 

In  this  work,  we  addressed  a  much  more  general  problem  than  what  has 
previously  been  addressed  in  the  literature.  Rather  than  limiting  the  study  to  a  simple 
problem  that  could  be  reduced  to  a  second-order  system  or  single-input  control  problem, 
in  this  work  we  investigated  a  complex  DC  electrical  distribution  system.  The  electrical 
distribution  system  studied  has  multiple  inputs.  Some  of  those  inputs  were  voltage 
sources,  while  others  were  current  sources.  Unlike  previous  studies,  the  system  under 
investigation  could  not  reasonably  be  approximated  by  a  second-order  system. 

A.  SIGNIFICANT  CONTRIBUTIONS 

1.  LQR-SM  Controller 

By  expanding  the  scope  of  the  study  to  a  complex  MVDC  distribution  system,  we 
were  forced  to  employ  a  multi-input  control  scheme.  The  control  scheme  developed  is  an 
adaptive,  multi-rate  LQR  controller  that  uses  selected  cost  function  matrices.  The 
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proposed  solution  reduces  computational  complexity  and  improved  robustness  as 
compared  to  periodic-discrete  multi-rate  LQR. 

The  proposed  LQR-SM  control  scheme  described  in  this  work  was  demonstrated 
through  MATLAB  simulation  to  perform  with  near  equivalence  to  a  traditional  periodic- 
discrete  implementation  of  multi-rate  LQR.  The  LQR-SM  method  of  multi-rate  LQR 
produced  nearly  identical  control  inputs  and  system  regulation  as  compared  to  the 
classical  periodic-discrete  implementation  but  did  so  using  a  small  fraction  of  the 
computing  resources.  This  was  done  by  rotating  through  a  set  of  either  B  or  R  matrices  to 
solve  for  sub-cycle  feedback  gains  for  each  possible  combination  of  inputs  rather  than 
constructing  large  block-cyclic  matrices.  By  reducing  the  size  of  the  matrices  operated  on 
by  the  Riccati  equation  solver,  we  reduced  the  computational  load  drastically.  This 
reduction  on  computation  load  enabled  the  control  designer  to  employ  an  adaptive  LQR 
controller  rather  than  be  constrained  to  off-line  computation  only. 

A  further  advantage  of  the  LQR-SM  control  scheme  was  improved  scalability  and 
reliability.  In  Chapter  IV  Section  D,  we  considered  a  traditional  multi-rate  LQR 
implemented  on  a  thirteenth-order  system.  Using  a  thirteenth-order  system  with  an  eight- 
step  block-cycle  produced  enonnous  block-cyclic  A,  B,  Q,  and  R  matrices.  The 
MATLAB  Riccati  solver  dareQ  issued  several  warnings  regarding  ill-conditioned 
matrices.  The  dareQ  solver  was  pushed  past  its  limit  when  attempting  to  solve  for  a 
larger,  twentieth-order  circuit  model.  In  contrast,  LQR-SM  was  easily  scaled  to  a 
twentieth-order  circuit  system  with  16  steps  per  cycle  without  issuing  any  warnings. 

While  LQR-SM  was  investigated  for  a  naval  MVDC  distribution  system,  the 
technique  described  is  not  constrained  to  one  specific  class  of  problems.  The  LQR-SM 
controller  can  be  applied  broadly.  Whether  the  system  under  study  is  single  input  or 
multi-input,  small  or  complex,  fully  controllable  or  not,  LQR-SM  may  be  applied.  The 
only  restrictions  are  that  the  system  is  stabilizable,  Q  is  positive  semi-definite,  and  R  is 
positive  definite. 
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2.  Design  Process 

The  next  contribution  is  an  easily  understood  design  process.  While  there  are 
rules  of  thumb  regarding  LQR  controller  design,  they  are  inadequate  for  high-order 
problems.  When  dealing  with  large  and  complex  systems,  the  interaction  between  state- 
variables  and  the  LQR  Q  and  R  matrix  values  can  become  overwhelming. 

In  our  design  process,  the  circuit  specifications  and  controller  parameters  were 
determined  simultaneously.  Starting  with  a  passively  stable  distribution  system,  the 
genetic  algorithm  iterates  through  the  chosen  number  of  generations  to  select  circuit  and 
controller  parameter  combinations  which  optimize  a  fitness  function.  The  genetic 
algorithm  design  process  described  in  this  dissertation  utilized  time-domain  simulations 
of  a  worst-case  power  transient.  The  simulation  trial  results  were  measured  and  checked 
for  compliance  with  overshoot,  undershoot,  settling  time,  and  hannonic  content 
specifications.  Finally,  system  stored  energy  was  measured,  with  the  optimal  candidate 
having  the  lowest  energy.  Other  engineers  with  different  or  more  complex  objectives  may 
include  those  measurements  into  their  own  fitness  functions. 

Due  to  the  CPL  in  our  system,  it  was  imperative  that  we  begin  with  a  passively 
stable  system.  Because  the  system  we  studied  was  not  globally  stable,  we  could  not 
choose  random  values  for  the  genomes  of  the  initial  generation  of  our  genetic  algorithm. 
Consequently,  random  selection  or  even  poor  selection  of  the  initial  generation  of 
genomes  may  result  in  an  optimization  that  never  converges  to  a  solution  within  the 
ROA. 


3.  Performance  Trade-Offs 

Lastly,  we  explored  perfonnance  trade-offs.  Our  system  had  many  state  variables 
that  closely  tracked  with  each  other  because  impedances  were  small.  We  experimented 
with  reduced-order  controllers  and  with  a  non-adaptive  controller  to  detennine  what 
trade-offs,  if  any,  were  experienced  by  reducing  the  fidelity  of  the  controller  or  by 
considering  a  non-adaptive  approach.  For  this  study,  we  learned  that  reduced  fidelity 
controllers  that  ignored  small  impedances  had  reduced  computational  load  without  any 
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penalty  in  minimum  energy  storage  or  transient  dynamic  performance.  Even  the  non- 
adaptive  controller  perfonned  as  well  as  the  adaptive  controllers. 

The  major  differences  between  the  controllers  were  in  the  sizes  of  their  respective 
ROAs.  In  the  minimum  energy  configuration,  ROAs  varied  quite  widely,  even  though  the 
zone  buck  capacitors  used  in  each  controller  scenario  were  similar.  An  experiment  where 
controllers  used  identical  circuit  parameters  and  Q  and  R  matrix  values  revealed  that 
choice  of  controller  Q  and  R  values  plays  a  significant  role  in  the  size  of  ROAs. 

B.  FUTURE  WORK 

So  far,  LQR-SM  controllers  have  only  been  implemented  in  MATLAB  software 
for  an  average  value  model.  The  robustness  and  utility  of  the  control  method  could  be 
further  proved  through  implementation  on  a  test  bed.  Testing  the  control  scheme  on  a  test 
bed,  first  as  a  switch-mode  model  and  then  in  hardware,  could  prove  the  real-world 
applicability  of  this  method.  Furthermore,  application  to  a  physical  model  would  also 
allow  for  testing  the  controller’s  sensitivity  to  microcontroller  delay  time  and  sensor 
measurement  noise.  For  this  to  occur,  a  HESS  controlled  current  source  needs  to  be 
developed. 

In  addition,  an  investigation  into  ROA  size  is  warranted.  Since  a  system  this 
complex  defies  analytical  analysis  of  the  ROA,  a  computational  investigation  into  the 
factors  which  most  strongly  influence  ROA  size  would  be  instructive.  The  genetic 
algorithm  could  be  modified  to  assess  ROA  size  and  an  ROA  score  included  in  the  fitness 
function. 
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APPENDIX  A.  MATLAB  CODE 


1.  Simulated  Annealing 

%Synthetic  Annealing  Algorithm 

%This  algorithm  randomly  modifies  the  "best  performing"  iteration  of 
system 

%variables.  Initially  the  variance  of  change  to  variables  is  large,  but 

%then  gradually  reduces  to  settle  on  a  final  value. 

clear; 

Best_Score  =  le6;  %initial  best  score.  Lower  score  is  better. 

To  =  100;  %initial  system  variance 
T(l)  =  To;  %System  temp 
k  =  1; %number  of  iterations 
iter  =  1; 

run  Circuit_Data25; 

C  =  [Cbus  Cdl  Cd2  Czl  Cz2  Cz3  Cd3  Cbl  5*Cb2  Cb3]; 

Q  =  [11111111111111111111];  %Q  diagonals 
R  =  [1  1  le3  1  1];  %Generator  in  use  penalty,  input  unavailable 
penalty,  Isl  penalty,  Is2  penalty 

C_best  =  C;  Q_best  =  Q;  R_best  =  R; 

L  =  [Lgl  Lg2  Lg3  Lg4  Lzl  Lz2  Lz3  Lbl  Lb2  Lb3]; 
r  =  [Rgl  Rg2  Rg3  Rg4  Rzl  Rz2  Rz3  Rdl  Rd2  Rd3] ; 

V  =  [Vref  Vrefl  Vref2  Vref3]; 

minValC  =  [le-6  le-7  le-7  le-7  Czl  Cz2  Cz3  Cripplel  Cripple2  Cripple3]; 
maxValC  =  C; 
minValQ  =  Q; 
maxValQ  =  2A11*Q; 
minValR  =[11  le3  0.05  0.05]; 
maxValR  =  [2A4  2A4  le3  1  1]; 
minVal  =  [minValC  minValQ  minValR] ; 
maxVal  =  [maxValC  maxValQ  maxValR] ; 

J(l)  =  fitness (V,  r,  L,  [C_best  Q_best  R_best]); 

J_best  =  J ( 1 ) ; 

while  T  >  0 
k  =  k+1 ; 

%Cooling  Profile 
if  k  >  2800 

T (k)  =  T(k-l)  -  40/400; 
elseif  k  >  1600 
T(k)  =  T(k-l)  -  60/1200; 
else 

T (k)  =  T(k-l) ; 
end 

%Available  variables  for  modification  are  bus  capacitance, 
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zone  input 


%filter  capacitances,  buck  capacitances,  R-matrix  diagonal  values, 

%  and  Q-matrix  diagonal  values. 

neighbor  =  mutateSA ( [C_best  Q_best  R_best] ,  minVal,  maxVal,  T (k)  ,  To); 

J_neighbor  =  fitness (V,  r,  L,  neighbor) ; 

if  J_neighbor  <  J_best 
J_best  =  J_neighbor; 

C_best  =  neighbor ( 1 : 10 ) ; 

Q_best  =  neighbor ( 11 : 30 ) ; 

R_best  =  neighbor ( 31 : 35 ) ; 
end 

J(k)  =  J_best; 

[k  T  (k)  J_best] 
end 


function  [  y  ]  =  mutateSA (  seed,  minVal,  maxVal,  T,  To  ) 

%MUTATE  returns  a  scalar  value  based  on  a  random  normal  distribution 
with  mean  of  'seed',  variance  of  0.1  and  with  cut-off  minimum  and 
maximum  values. 

%This  helps  to  choose  only  20%  of  genes  for  mutation 
chooser  =  5*rand(l,  length ( seed) ) ; 

for  k  =  1 : length ( seed) 
if  chooser  (k)  >  1 

y (k)  =  seed (k) ; 
else 

y(k)  =  min (max (minVal (k) ,  seed(k)*(l  +  0 . 5*T/To*randn (1) )  ), 

maxVal (k) ) ; 

end 

end 

%enforce  that  small  generators  and  large  generators  have  same  penalties 
y (13)  =  y  (11) ; 
y (14)  =  y (12)  ; 
end 


2.  Genetic  Algorithm 

3.  Super-Routine 

%%%Genetic  Algorithm  -  GA  with  very  tight  restrictions  on  C  Q  R%%% 

%This  script  inputs  a  circuit  parameters  file,  including  all  of  the  RLC 
%data  for  a  Four-Generator,  Three-Zone  MVDC  shipboard  distribution 
system 

%with  a  common  MVDC  bus. 
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%Since  the  system  will  be  controlled  via  an  adaptively  linearizing  LQR 
%controller,  the  Q  and  R  matrices  will  be  inputs  as  well. 

%Start  with  a  system  that  is  passively  stable,  then  selectively  reduce 
%capacitor  sizes,  increase  variable  error  penalties,  and  reduce  input 
%penalties  to  produce  a  system  with  acceptable  performance  and  minimum 
%stored  energy. 

clear; 

%clc; 

%Import  Circuit  Data  for  the  passively  stable  system 
run  Circuit_Data25; 

C  =  [2*Cbus  Cdl  Cd2  Cd3  Czl  Cz2  Cz3  Cbl  5*Cb2  Cb3]; 

L  =  [Lgl  Lg2  Lg3  Lg4  Lzl  Lz2  Lz3  Lbl  Lb2  Lb3]; 
r  =  [Rgl  Rg2  Rg3  Rg4  Rzl  Rz2  Rz3  Rdl  Rd2  Rd3]  ; 

V  =  [Vref  Vrefl  Vref2  Vref3]; 

ReL  =  1;  %Input  penalty  for  large  generators. 

ReS  =  1;  %Input  penalty  for  small  generators. 

Pen  =  le3;  %Penalty  value  for  unusable  inputs. 

Rsl  =  l;%Isl  penalty 
Rs2  =  1;%Is2  penalty 

%Variable  Error  Penalty  Matrix 
Q  =  ones (1,20) ; 

%Input  Penalty  Matrix 

R  =  [ReL  ReS  Pen  Rsl  Rs2];  %Input  penalty  for  large  generator,  small 
generator,  unavailable  input,  Isl,  Is2 

%Penalty  values 
os_penalty  =  50; 
fft_penalty  =  50; 

%Define  the  seed/best 
C_best  =  C; 

Q_best  =  Q; 

R_best  =  R; 

J_best  =  101; 

%Create  Generation  #1 
pop  =  16;  %population 
for  i  =  1 : pop 
if  i  ==  1 

G  (i,  :  )  =  [C  Q  R]  ; 

else 

minValC  =  [le-6  le-7  le-7  le-7  Czl  Cz2  Cz3  Cripplel  Cripple2  Cripple3] ; 
maxValC  =  C; 
minValQ  =  Q; 
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maxValQ  =  2~11*Q; 
minValR  =[11  le3  0.05  0.05]; 
maxValR  =  [2A4  2A4  le3  1  1]; 
minVal  =  [minValC  minValQ  minValR] ; 
maxVal  =  [maxValC  maxValQ  maxValR] ; 

G(i,  :)  =  mutateO ( [C  Q  R] ,  minVal,  maxVal); 

end 

end 

Generations  =  200; 
for  k  =  1 : Generations 
%tic 

%Test  fitness  of  each  member  of  the  generation 
for  i  =  1 : pop 

Jgen(i)  =  fitness (V,  r,  L,  G(i,  :)  ); 

end 

J (k)  =  min ( Jgen) 

%Determine  which  population  members  are  breeders.  Here  we  use  top  half 
m  =  1 ; 

for  i  =  1 : pop 
if  Jgen(i)<  median (Jgen) 

Breeder (m,  :)  =  G(i,  :); 

Jbreed(m)  =  Jgen(i); 

m  =  m+1; 

end 

end 

%Use  parents  to  create  children.  Parents  survive  to  the  next 
%generation 

for  i  =  1 : length (Jbreed) /2 

mixerl  =  round (rand (1,  length ( [C  Q  R] ) ) ) ; 

mixer2  =  1 -mixerl; 

parentl  =  Breeder (2*i-l,  :); 

parent2  =  Breeder (2*i,  :); 

childl  =  parentl . *mixerl  +  parent2 . *mixer2 ; 
child2  =  parentl . *mixer2  +  parent2 . *mixerl ; 
if  Jbreed (2*i-l )  ==  min (Jgen) 

%do  nothing  -  this  is  the  elite  parent 
C_best  =  parentl (1 : 10) ; 

Q_best  =  parentl ( 11 : 30 ) ; 

R_best  =  parentl ( 31 : 35 ) ; 
else 

%mutate  parentl 

parentl  =  mutate (parentl ,  minVal,  maxVal); 
end 

if  Jbreed (2*i)  ==  min (Jgen) 

%do  nothing  -  this  is  the  elite  parent 
C_best  =  parent2 ( 1 : 10 ) ; 

Q_best  =  parent2 (11 : 30) ; 

R_best  =  parent2 ( 31 : 35 ) ; 
else 

%mutate  parent2 

parent2  =  mutate (parent2 ,  minVal,  maxVal); 
end 
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childl  =  mutate (childl,  minVal,  maxVal); 
child2  =  mutate (child2,  minVal,  maxVal); 

G(2*i-1,  :)  =  parentl; 

G(2*i,  :)  =  parent2; 

G(pop/2  +  2*i-l,  :)  =  childl; 

G(pop/2  +  2*i,  :)  =  child2; 

end 

%toc 

end 


4.  Mutate  Function 

function  [  y  ]  =  mutateO (  seed,  minVal,  maxVal  ) 

%MUTATE  returns  a  scalar  value  based  on  a  random  normal  distribution 
with  mean  of  'seed',  variance  of  0.1  and  with  cut-off  minimum  and 
maximum  values, 
for  k  =  1 : length ( seed) 

y(k)  =  minVal (k)  +  (maxVal (k)  -  minVal (k) ) *rand (1) ;  %min (max (minVal (k) , 

seed(k)*(l  +  0 . 5*randn ( 1 ) )  ),  maxVal (k) ) ; 

end 

%enforce  that  small  generators  and  large  generators  have  same  penalties 
y (13)  =  y (11) ; 
y (14)  =  y (12)  ; 

end 


function  [  y  ]  =  mutate (  seed,  minVal,  maxVal  ) 

%MUTATE  returns  a  scalar  value  based  on  a  random  normal  distribution 
with  mean  of  'seed',  variance  of  0.1  and  with  cut-off  minimum  and 
maximum  values. 

%This  helps  to  choose  only  20%  of  genes  for  mutation 
chooser  =  5*rand(l,  length ( seed) ) ; 

for  k  =  1 : length ( seed) 
if  chooser (k)  >  1 

y (k)  =  seed (k) ; 
else 

y(k)  =  min (max (minVal (k) ,  seed(k)*(l  +  0 . 5*randn ( 1 ) )  ),  maxVal (k) ) ; 

end 

end 

%enforce  that  small  generators  and  large  generators  have  same  penalties 
y (13)  =  y (11) ; 
y (14)  =  y (12)  ; 
end 
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5, 


Fitness  Function 


function  [  J__trial  ]  =  fitness  (  V,  r,  L,  G  ) 

%Run  a  trial  on  each  population  member,  then  return  teh  fitness  score 

%Detailed  explanation  goes  here 

%Run  a  Trial 

run  Circuit_Data25 

os_penalty  =  50; 
fft_penalty  =  50; 

C  =  G  (1 : 10)  ; 

Q  =  G  (11 :  30)  ; 

R  =  G ( 31 : 35 ) ; 

H  =  trial20 (V,  r,  L,  C,  Q,  R)  ; 

Vbus  =  H (1,  :  ) ; 

Vzl  =  H  (2,  :  )  ; 

Vz2  =  H (3,  :  )  ; 

Vz3  =  H  (4,  :  )  ; 

Vbl  =  H  (5,  :  )  ; 

Vb2  =  H  (6,  :  )  ; 

Vb3  =  H  (7,  :  )  ; 

El  =  H  (8,  :  )  ; 

E2  =  H  (9,  :  )  ; 

E3  =  H  (10,  :  )  ; 

E4  =  H  ( 1 1 ,  :)  ; 

Isl  =  H  (12,  :  )  ; 

Is2  =  H  (13,  :  )  ; 

%%Score  the  trial 


%Overshoot  score 
J_overshoot  =  0; 


overshoot_Vbus  =  max (  (max (Vbus) -Vref) 
%(  max (Vbus)  -  min (Vbus)  )/Vref; 
overshoot_Vbl  =  max (  (max (Vbl ) -Vref 1 ) , 
%(  max (Vbl)  -  min (Vbl)  )/Vrefl; 
overshoot_Vb2  =  max (  (max (Vb2 ) -Vref 2 ) , 
%(  max(Vb2)  -  min(Vb2)  ) /Vref 2; 
overshoot_Vb3  =  max (  (max (Vb3 ) -Vref 3 ) , 
%(  max(Vb3)  -  min(Vb3)  ) /Vref3; 


abs (min (Vbus ) 
abs (min (Vbl )  - 
abs (min (Vb2 )  - 

abs (min (Vb3 )  - 


-  Vref)  ) /Vref ; 
Vref 1 )  ) /Vref 1 ; 
Vref2)  )/Vref2; 
Vref 3 )  ) /Vref 3 ; 


if  overshoot_Vbus  >  0.10 

J_overshoot  =  J_overshoot  +  os_penalty; 

end 

if  overshoot_Vbl  >  0.10 

J_overshoot  =  J_overshoot  +  os_penalty; 
end 

if  overshoot_Vb2  >  0.10 

J_overshoot  =  J_overshoot  +  os_penalty; 
end 

if  overshoot_Vb3  >  0.10 
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J_overshoot  =  J_overshoot  +  os_penalty; 
end 


%Harmonic  content  score 
J_harmonics  =  0; 
if  sldetect (Vbus )  >  -60 
J_harmonics  =  J_harmonics  + 
end 

if  sldetect (Vbl )  >  -60 
J_harmonics  =  J_harmonics  + 
end 

if  sldetect (Vb2 )  >  -60 
J_harmonics  =  J_harmonics  + 
end 

if  sldetect (Vb3 )  >  -60 
J_harmonics  =  J_harmonics  + 
end 


f ft_penalty ; 


f ft_penalty ; 


f ft_penalty ; 


f ft_penalty ; 


%Capacitance  Score 

J_capacitance  =  sum((C  -  [0000  Czl  Cz2  Cz3  Cripplel  Cripple2 
Cripple3] )  . * [1  1  1  1  0  0  0  (Vrefl/Vref)  (Vref2/Vref)  (Vref3/Vref) ] )  ; 


%Stored  Energy 

E_caps  =  (1/2) *C* [Vref  Vref  Vref  Vref  Vref  Vref  Vref  Vrefl  Vref2 
Vref 3 ] . A2 ' ; 


E_Isl  =  zeros (size (Isl)  )  ; 

E_Is2  =  zeros (size (Is2) ) ; 
dt  =  0 . 100/ (length (Isl)  -  1); 

for  k  =  1 : length ( Isl ) 

E_Isl (k)  =  sum (  Vbl (1 :k) . *Isl (1 :k)  )*dt; 

E_Is2 (k)  =  sum (  Vb2 (1 : k) . *Is2 (1 : k)  )*dt; 
end 

E_Isl_max  =  max(E_Isl); 

E_Is2_max  =  max(E_Is2); 

J_Energy  =  (E_caps  +  E_Isl_max  +  E_Is2_max) /le6; 

%Trial  Failure  Score 
J_t  f  =  0 ; 

if  length (Vbus)  <  4000 

J_tf  =  100; 

end 


%J_summary  =  [J_overshoot  J_harmonics  J_tf  J_Energy] 

J_trial  =  J_overshoot  +  J_harmonics  +  J_tf  +  J_Energy;  % J_capacitance ; 


end 


109 


6, 


Trial  Function 


funct. 

ion 

[YY] 

=  tri 

.al20 (Vi,  ri. 

Li 

,  Ci, 

Qi, 

Ri) 

Vref  : 

=  Vi  (1)  ; 

Vrefl 

Vi 

(2);  Vref 2 

=  Vi 

(3)  ; 

Vref 3  =  Vi ( 4 ) ; 

Rgl  = 

ri 

(1)  ; 

Rg2  = 

ri 

(2)  ; 

:  Rg3  = 

ri 

(3)  ; 

Rg4  = 

=  ri  (4) ; 

Rz  1  = 

ri 

(5)  ; 

Rz2  = 

ri 

(6)  ; 

:  Rz3  = 

ri 

(7)  ; 

Rdl  = 

ri 

(8)  ; 

Rd2  = 

ri 

(9)  i 

:  Rd3  = 

ri 

(10)  ; 

Lgl  = 

Li 

(1)  ; 

Lg2  = 

Li 

(2)  ; 

:  Lg3  = 

Li 

(3)  ; 

Lg4  = 

=  Li  (4)  ; 

Lzl  = 

Li 

(5)  ; 

Lz2  = 

Li 

(6)  ; 

:  Lz3  = 

Li 

(7)  ; 

Lbl  = 

Li 

(8)  ; 

Lb  2  = 

Li 

(9)  i 

:  Lb  3  = 

Li 

(10)  ; 

Cbus  : 

=  C 

i  (1)  ; 

Cdl  = 

=  C 

i  (2) 

i  ;  Cd2  = 

=  C 

i  (3)  ; 

Cd3 

=  Ci  (4)  ; 

Czl  = 

Ci 

(5)  ; 

Cz2  = 

Ci 

(6)  ; 

:  Cz3  = 

Ci 

(7)  ; 

Cbl  = 

Ci 

(8)  ; 

Cb2  = 

Ci 

(9)  i 

:  Cb3  = 

Ci 

(10)  ; 

Tesd  =  1/I6e3; 

Tgen  =  l/4e3; 

%State-Variable  Initial  values 


%Circuit_Data25 

PI  =  [15  5  28]*le6;  %Initial  Power  for  zones  1,2,3 
P2  =  [20  30  46]*le6;  %Final  Power 

P  =  PI; 

R  =  [Vrefl  Vref2  Vref3] . A2./P; 

%Generator  Currents 

Iglkml  =  0 . 40*sum (P) /Vref ; %Generator  current 
Ig2kml  =  0 . 10*sum (P) /Vref ; %Generator  current 
Ig3kml  =  0 . 40*sum (P) /Vref ; %Generator  current 
Ig4kml  =  0 . 10*sum (P) /Vref ; %Generator  current 
Itotalkml  =  Iglkml+Ig2kml+Ig3kml+Ig4kml ; 


%Zone  Down-stepped  values 
Iblkml  =  P(l) /Vrefl; 
Vblkml  =  Vrefl; 

Vblkm2  =  Vblkml; 


Ib2kml  =  P(2)/Vref2; 
Vb2kml  =  Vref 2; 
Vb2km2  =  Vb2kml; 


Ib3kml  =  P(3)/Vref3; 
Vb3kml  =  Vref 3; 
Vb3km2  =  Vb3kml; 


%Cable  Currents 

Izlkml  =  (Vref 1/Vref ) *Iblkml; %Load  current  zone  1 
Iz2kml  =  (Vref2/Vref ) *Ib2kml; %Load  current  zone  2 
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Iz3kml 


(Vref 3/Vref ) *Ib3kml; %Load  current  zone  3 


%Bus  Voltage 

Vbuskm2  =  Vref;%MVDC  bus  voltage 
Vbuskml  =  Vbuskm2; 


%Point 

of  Load  Voltages 

Vzlkm2 

= 

Vref ; %Voltage  at  input 

to  zone#l 

CPL 

Vzlkml 

= 

Vzlkm2 ; 

Vdlkm2 

= 

Vref;%Damper  capacitor 

voltage 

Vdlkml 

= 

Vdlkm2 ; 

Vz2km2 

= 

Vref ; %Voltage  at  input 

to  zone#2 

CPL 

Vz2kml 

= 

Vz2km2 ; 

Vd2km2 

= 

Vref;%Damper  capacitor 

voltage 

Vd2kml 

= 

Vd2km2 ; 

Vz3km2 

= 

Vref ; %Voltage  at  input 

to  zone#3 

CPL 

Vz3kml 

= 

Vz3km2 ; 

Vd3km2 

= 

Vref;%Damper  capacitor 

voltage 

Vd3kml 

= 

Vd3km2 ; 

%Input 

Device  Values 

Elkm2  =  Vref  t  Rgl*Iglkml; %Generator  voltage 
Elkml  =  Elkm2 ; 

E2km2  =  Vref  +  Rg2*Ig2kml; %Generator  voltage 
E2kml  =  E2km2; 

E3km2  =  Vref  +  Rg3*Ig3kml ; %Generator  voltage 
E3kml  =  E3km2; 

E4km2  =  Vref  +  Rg4*Ig4kml; %Generator  voltage 
E4kml  =  E4km2; 

Islkml  =  0;%Zone  1 
Is2kml  =  0;%Zone  2 
Is3kml  =  0; 

dt  =  5e-9; 
t  =  0:dt:2*le7*dt; 

%decimate  =  length (t) /400; 

tau  =  0.05; 

Pzlestkml  =  P(l); 

Pz2estkml  =  P(2); 

Pz3estkml  =  P(3); 

Pestkml  =  sum(P); 

Rbl  =  1;  Rb2  =  1;  Rb3  =  1; 

dl  =  Vrefl/Vref;  d2  =  Vref2/Vref;  d3  =  Vref3/Vref; 


%Pen  =  le3;%Penalty  for  unusable  inputs; 
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index  =  1; 
decimate  =  5e3; 
token  =  0; 


for 

k  =  1 :  length (t ) 

Pzl 

=  Vblkml* 

(  Iblkml  -  Cbl* 

(Vblkml 

-  Vb 1 km2 ) / dt 

+ 

Islkml ) ; 

Pz2 

=  Vb2kml* 

(  Ib2kml  -  Cb2* 

(Vb2kml 

-  Vb2km2 ) / dt 

+ 

Is2kml ) ; 

Pz3 

=  Vb3kml* 

(  Ib3kml  -  Cb3* 

(Vb3kml 

-  Vb3km2 ) / dt 

+ 

Is3kml ) ; 

if  : 

mod (t (k) , 

Tesd)  ==  0 

%Implement  LQR 

Pzlestk  =  tau*Pzl  +  (1-tau) *Pzlestkml; %Zone  power  estimates 
Pz2estk  =  tau*Pz2  +  (1-tau) *Pz2estkml; 

Pz3estk  =  10*tau*Pz3  +  (l-10*tau) *Pz3estkml; 

Pestk  =  Pzlestk  +  Pz2estk  +  Pz3estk; 

10  =  Pestk/Vref;  %Equilibrium  current  value 

101  =  Pzlestk/Vref ; 

102  =  Pz2estk/Vref ; 

103  =  Pz3estk/Vref ; 

Iobl  =  Pz lestk/Vref 1 ; 

Iob2  =  Pz2estk/Vref 2 ; 

Iob3  =  Pz3estk/Vref 3; 

%CPL  non-linear  resistance  values 

Rnll  =  -VblkmlA2/Pzlestk;  %-Vzlkml A2/Pzlestk; 

Rnl2  =  -Vb2kmlA2/Pz2estk;  %-Vz2kml A2 /Pz2estk; 

Rnl3  =  -Vb3kmlA2/Pz3estk;  %-Vz3kml A2/Pz3estk; 


XI  = 

Iglkml  - 

. 40*10; 

X2  = 

Ig2kml  - 

.10*10; 

X3  = 

Ig3kml  - 

. 40*10; 

X4  = 

Ig4kml  - 

.10*10; 

X5  = 

Vbuskml 

-  Vref; 

X6  = 

Izlkml  - 

Iol; 

X7  = 

Iz2kml  - 

Io2 ; 

X8  = 

Iz3kml  - 

Io3  ; 

X9  = 

Vzlkml  - 

Vref  - 

X10 

=  Vz2kml 

-  Vref  - 

Xll 

=  Vz3kml 

-  Vref  - 

X12 

=  Vdlkml 

-  Vref  - 

X13 

=  Vd2kml 

-  Vref  - 

XI 4 

=  Vd3kml 

-  Vref  - 

X15 

=  Iblkml 

-  Iobl; 

XI 6 

=  Ib2kml 

-  Iob2; 

X17 

=  Ib3kml 

-  Iob3; 

X18 

=  Vblkml 

-  Vrefl; 

XI 9 

=  Vb2kml 

-  Vref2; 

X20 

=  Vb3kml 

-  Vref3; 

%Change  of  variables 


Rzl*Iol; 
Rz2*Io2 ; 
Rz3*Io3; 
Rzl*Iol ; 
Rz2*Io2 ; 
Rz3*Io3; 
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X  =  [XI;  X2 ;  X3;  X4 ;  X5;  X6;  X7;  X8 ;  X9;  X10;  Xll;  X12;  X13;  X14 
XI 6 ;  X17 ;  X18;  X19;  X20]; 


%ESDs  switching  but  generators  are  not 
if  mod(t(k),  Tgen)  ~=  0 
mode  =  0 ; 

U  =  SysControl20 (X, dl, d2, d3,  Rnll, Rnl2, Rnl3,  Li,  ri,  Ci,  Qi,  Ri, 
Islk  =  U  (18) ; 

Is2k  =  U  (19) ; 

Is3k  =  U  (20) ; 
end 

%  ESDs  and  One  generator  are  switching  together 
if  mod(t(k),  Tgen)  ==  0 
if  token  ==  0 
mode  =  1 ; 

U  =  SysControl20 (X, dl, d2, d3,  Rnll, Rnl2, Rnl3,  Li,  ri,  Ci,  Qi,  Ri, 
Elk  =  U(l)  +  Vref  +  Rgl*0. 40*10; 

E2k  =  E2kml;  %U(2)  +  Vref  t  Rg2*0. 10*10;  %E2kml; 

E3k  =  E3kml ;  %U(3)  +  Vref  +  Rg3*0. 40*10;  %E3kml; 

E4k  =  E4kml ;  %U(4)  +  Vref  t  Rg4*0. 10*10;  %E4kml; 

Islk  =  U  (18) ; 

Is2k  =  U  (19) ; 

Is3k  =  U  (20) ; 
token  =  1; 
elseif  token  ==  1 
mode  =  2 ; 

U  =  SysControl20 (X, dl, d2, d3,  Rnll, Rnl2, Rnl3,  Li,  ri,  Ci,  Qi,  Ri, 
Elk  =  Elkml ;  %U(1)  +  Vref  t  rgl*0 . 40*10 (k) ; 

E2k  =  U (2 )  +  Vref  +  Rg2*0. 10*10; 

E3k  =  E3kml;  %U(3)  +  Vref  +  rg3*0 . 40*10 (k) ; 

E4k  =  E4kml ;  %U(4)  +  Vref  +  rg4*0 . 10*10 (k) ; 

Islk  =  U (18) ; 

Is2k  =  U (19) ; 

Is3k  =  U (20) ; 
token  =  2; 
elseif  token  ==  2 
mode  =  3; 

U  =  SysControl20  (X,  dl,  d2,  d3,  Rnll,  Rnl2,  Rnl3,  Li,  ri,  Ci,  Qi,  Ri, 
Elk  =  Elkml;  %U(1)  +  Vref  +  rgl*0 . 40*10 (k) ; 

E2k  =  E2kml;  %U(2)  +  Vref  t  rg2*0 . 10*10 (k) ; 

E3k  =  U ( 3 )  +  Vref  +  Rg3*0. 40*10; 

E4k  =  E4kml ;  %U(4)  +  Vref  t  rg4*0 . 10*10 (k)  ; 

Islk  =  U  (18) ; 

Is2k  =  U  (19) ; 

Is3k  =  U  (20) ; 
token  =  3; 
elseif  token  ==  3 
mode  =  4 ; 

U  =  SysControl20 (X, dl, d2, d3,  Rnll, Rnl2, Rnl3,  Li,  ri,  Ci,  Qi,  Ri, 
Elk  =  Elkml;  %U(1)  +  Vref  t  rgl*0 . 40*10 (k) ; 

E2k  =  E2kml;  %U(2)  +  Vref  +  rg2*0 . 10*10 (k) ; 

E3k  =  E3kml;  %U(3)  +  Vref  t  rg3*0 . 40*10 (k) ; 

E4k  =  U ( 4 )  +  Vref  +  Rg4*0. 10*10; 

Islk  =  U  (18) ; 


;  XI 5; 


mode ) 


mode ) 


mode ) 


mode ) 


mode ) 
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Is2k  =  U (19) ; 
Is3k  =  U (20) ; 
token  =  0; 
end 
end 


end 


Is3k  =  0; 

%State-Space  equations 
Idlk  =  (Vzlkml  -  Vdlkml)/Rdl; 
Id2k  =  (Vz2kml  -  Vd2kml)/Rd2; 
Id3k  =  (Vz3kml  -  Vd3kml)/Rd3; 


dlgl  = 
dlg2  = 
dlg3  = 
dlg4  = 
dVbus  = 
Iz3kml ) 
dlzl  = 
dlz2  = 
dlz3  = 
dVzl  = 
dVz2  = 
dVz3  = 
dVdl  = 
dVd2  = 
dVd3  = 
dlbl  = 
dlb2  = 
dlb3  = 
dVbl  = 
dVb2  = 
dVb3  = 


(l/Lgl)*(Elk  -  Rgl*Iglkml  -  Vbuskml); 

(l/Lg2)*(E2k  -  Rg2*Ig2kml  -  Vbuskml); 

(l/Lg3)*(E3k  -  Rg3*Ig3kml  -  Vbuskml); 

(l/Lg4)*(E4k  -  Rg4*Ig4kml  -  Vbuskml); 

(1/Cbus) * (Iglkml  +  Ig2kml  +  Ig3kml  + 


Ig4kml  -  Izlkml  - 


(1/Lzl) 
(1/Lz2) 
(1/Lz3) 
(1/Czl) 
( 1/Cz2 ) 
(1/Cz3) 
( 1/Cdl ) 
(1/Cd2) 
(1/Cd3) 
(1/Lbl) 
(1/Lb2) 
(1/Lb3) 
(1/Cbl) 
( 1/Cb2 ) 
(1/Cb3) 


* (Vbuskml  - 
* (Vbuskml  - 
* (Vbuskml  - 
* (  Izlkml  - 
* (  Iz2kml  - 
* (  Iz3kml  - 
*Idlk; 

* Id2k; 
*Id3k; 

* (dl*Vzlkml 
* (d2*Vz2kml 
* (d3*Vz3kml 
* (Iblkml 
* ( Ib2kml 
* ( Ib3kml 


Rzl*Izlkml 
Rz2*Iz2kml 
Rz3*Iz3kml 
Idlk  -  dl* 
Id2k  -  d2* 
Id3k  -  d3* 


-  Vzlkml) 

-  Vz2kml) 

-  Vz3kml) 
Iblkml ) ; 
Ib2kml ) ; 
Ib3kml ) ; 


-  Vblkml ) ; 

-  Vb2kml ) ; 

-  Vb3kml ) ; 
P (1) /Vblkml 
P (2) /Vb2kml 
P (3) /Vb3kml 


+  Islk  ) 
+  Is2k  ) 
+  Is3k  ) 


%Variable  updates 

Iglk  =  Iglkml  +  dlgl*dt; 

if  Iglk  <  0 

Iglk  =  0; 

end 

Ig2k  =  Ig2kml  +  dlg2*dt; 
if  Ig2k  <  0 
Ig2k  =  0; 
end 

Ig3k  =  Ig3kml  +  dlg3*dt; 
if  Ig3k  <  0 
Ig3k  =  0; 
end 


Ig4k  =  Ig4kml  +  dlg4*dt; 
if  Ig4k  <  0 
Ig4k  =  0; 


Iz2kml  - 
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end 

Vbusk 

=  Vbuskml 

+  dVbus*dt; 

Izlk 

= 

Izlkml 

+ 

dlzl*dt; 

Iz2k 

= 

Iz2kml 

+ 

dlz2*dt ; 

Iz3k 

= 

Iz3kml 

+ 

dlz3*dt; 

Vzlk 

= 

Vzlkml 

+ 

dV  z 1  *  dt ; 

Vz2k 

= 

Vz2kml 

+ 

dVz  2  *  dt ; 

Vz3k 

= 

Vz3kml 

+ 

dVz  3  *  dt ; 

Vdlk 

= 

Vdlkml 

+ 

dVdl*dt ; 

Vd2k 

= 

Vd2kml 

+ 

dVd2*dt; 

Vd3k 

= 

Vd3kml 

+ 

dVd3*dt ; 

Iblk 

= 

Iblkml 

+ 

dlbl*dt; 

if  Iblk  <  0 

Iblk 

= 

0; 

end 

Ib2k 

= 

Ib2kml 

+ 

dlb2*dt; 

if  Ib2k  <  0 

Ib2k 

= 

0; 

end 

Ib3k 

= 

Ib3kml 

+ 

dlb3*dt ; 

if  Ib3k  <  0 

Ib3k 

= 

0; 

end 

Vblk 

= 

Vblkml 

+ 

dVb  1  *  dt ; 

Vb2k 

= 

Vb2kml 

+ 

dVb2*dt; 

Vb3k 

= 

Vb3kml 

+ 

dVb  3  *  dt ; 

%Break  out  of  the  trial  if  it  appears  caret)  will  fail 

if  (Vblk/Vref 1 )  <0.1 

break 

elseif  (Vb2k/Vref2)  <  0.1 
break 

elseif  (Vb3k/Vref3)  <  0.1 

break 

end 


Vzlkm2  = 

Vzlkml 

Vz2km2  = 

Vz2kml 

Vz3km2  = 

Vz3kml 

Vblkm2  = 

Vblkml 

Vb2km2  = 

Vb2kml 

Vb3km2  = 

Vb3kml 

Iglkml  = 

Iglk; 

Ig2kml  = 

Ig2k; 

Ig3kml  = 

Ig3k; 

Ig4kml  = 

Ig4k; 

Izlkml  = 

Izlk; 

Iz2kml  = 

I  z2k; 

Iz3kml  = 

Iz3k; 

Vbuskml  = 

=  Vbusk 

Vzlkml  = 

Vzlk; 

Vz2kml  = 

Vz2k; 

Vz3kml  = 

Vz3k; 
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Vdlkml  =  Vdlk; 

Vd2kml  =  Vd2k; 

Vd3kml  =  Vd3k; 

Iblkml  =  Iblk; 

Ib2kml  =  Ib2k; 

Ib3kml  =  Ib3k; 

Vblkml  =  Vblk; 

Vb2kml  =  Vb2k; 

Vb3kml  =  Vb3k; 

Pzlestkml  = 

Pzlestk; 

Pz2estkml  = 

Pz2estk; 

Pz3estkml  = 

Pz3estk; 

Elkml  =  Elk 

r 

E2kml  =  E2k 

r 

E3kml  =  E3k 

r 

E4kml  =  E4k 

r 

Islkml  =  Is lk; 

Is2kml  =  Is2k; 

Is3kml  =  Is3k; 

%decimation 

routine 

if  mod(k,  decimate)  ==  0 

Igl ( index) 

=  Iglk; %Generator  Currents 

Ig2 ( index) 

=  Ig2k; 

Ig3 ( index) 

=  Ig3k; 

Ig4 ( index) 

=  Ig4k; 

Vbus ( index) 

=  Vbusk;%Main  bus  voltage 

Izl (index) 

=  Izlk;%Zone  cable  currents 

Iz2 (index) 

=  Iz2k; 

Iz3 (index) 

=  Iz3k; 

Vzl (index) 

=  Vz lk; %Voltages  at  input  to  buck  converter 

Vz2 ( index) 

=  Vz2k; 

Vz3 (index) 

=  Vz3k; 

Vdl ( index) 

=  Vdlk;%Damper  capacitor  voltages 

Vd2 ( index) 

=  Vd2k; 

Vd3 ( index) 

=  Vd3k; 

Idl ( index) 

=  Idlk;%Damper  capacitor  currents 

Id2 ( index) 

=  Id2k; 

Id3 ( index) 

=  Id3k; 

Ibl ( index) 

=  Iblk;%Buck  chopper  currents 

Ib2 ( index) 

=  Ib2k; 

Ib3 ( index) 

=  Ib3k; 

Vbl ( index) 

=  Vblk;%Buc  chopper  voltages 

Vb2 ( index) 

=  Vb2k; 

Vb3 ( index) 

=  Vb3k; 

El ( index)  = 

Elk; %Generator  VOltages 

E2 ( index)  = 

E2k; 

E3 ( index)  = 

E3k; 

E4 ( index)  = 

E4k ; 

Isl (index) 

=  Islk; %In jected  Currents 

Is2 (index) 

=  Is2k; 

Is3 (index) 

=  Is3k; 

Pest ( index) 

=  Pestk;%Total  Estimated  load  power 
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Pz lest ( index)  =  Pzlestkml ; %Zone  estimated  load  power 
Pz2est ( index)  =  Pz2estkml; 

Pz3est (index)  =  Pz3estkml; 

Iobl_ (index)  =  Iobl ; %Desired  buck  chopper  current 
Iob2_(index)  =  Iob2; 

Iob3_(index)  =  Iob3; 
time (index)  =  t  (k) ; 

%PP ( index)  =  P ; 

%Icpl (index)  =  P/Vzlk; 

%Vest (index)  =  Vestk; 
index  =  index  +  1; 
end 

%Change  load  power 

if  (  t (k)  >  0.001)&&(  t (k)  <  0.050  ) 

P  =  P2 ; 

R  =  [Vrefl  Vref2  Vref3] .A2./P; 
elseif  t (k)  >  0.050 
P  =  PI; 

R  =  [Vrefl  Vref2  Vref3] . A2./P; 
end 

end 

YY  =  [ Vbus ;  Vzl;  Vz2;  Vz3;  Vbl;  Vb2;  Vb3;  El;  E2;  E3;  E4;  Isl;  Is2]; 


7.  Controller  Function 


function  [U]  =  SysControl20 (X, dl, d2, d3,  Rnll, Rnl2, Rnl3,  Li,  ri. 


Ri, 

mode) ; 

%Th: 

Ls  functii 

on  to 

be 

used 

in 

conjuction  with 

CPL  case 

Rgl 

=  ri  (1)  ; 

Rg2 

=  ri 

(2)  ; 

Rg3 

=  ri ( 3 ) ;  Rg4  = 

ri  (4) ; 

Rzl 

=  ri  (5)  ; 

Rz2 

=  ri 

(6)  ; 

Rz3 

=  ri  (7) ; 

Rdl 

=  ri  (8)  ; 

Rd2 

=  ri 

(9)  ; 

Rd3 

=  ri  (10)  ; 

Ci,  Qi, 


Lgl  =  Li (1) ;  Lg2 
Lzl  =  Li  (5) ;  Lz2 
Lbl  =  Li  (8) ;  Lb2 


Li (2 ) ;  Lg3 
Li  (6) ;  Lz3 
Li  (9);  Lb 3 


Li  (3) ;  Lg4  =  Li  (4) ; 
Li  (7)  ; 

Li  (10)  ; 


Cbus  =  Ci (1)  ;  Cdl  =  Ci(2);  Cd2  =  Ci(3);  Cd3  =  Ci(4); 
Czl  =  Ci ( 5 ) ;  Cz2  =  Ci(6);  Cz3  =  Ci(7); 

Cbl  =  Ci  ( 8 ) ;  Cb2  =  Ci(9);  Cb3  =  Ci(10); 


Pen  =  Ri ( 3 ) ; 
if  mode  ==  0 
r_l  =  Pen; 
r_2  =  Pen; 
r_3  =  Pen; 
r_4  =  Pen; 

elseif  mode  == 
r_l  =  Ri  ( 1 )  ; 


1 
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r_2  =  Pen; 
r_3  =  Pen; 
r_4  =  Pen; 

elseif  mode  == 
r_l  =  Pen; 
r_2  =  Ri  (2  )  ; 
r_3  =  Pen; 
r_4  =  Pen; 

elseif  mode  ==  3 
r_l  =  Pen; 
r_2  =  Pen; 
r_3  =  Ri ( 1 ) ; 
r_4  =  Pen; 

elseif  mode  ==  4 
r__l  =  Pen; 
r_2  =  Pen; 
r_3  =  Pen; 
r_4  =  Ri  (2 )  ; 

end 

Bll  =  1;  B22  =  1;  B33  =  1;  B44  =  1; 

A  =  [-Rgl/Lgl  000  -1/Lgl  000000000000000;... 

0  -Rg2/Lg2  0  0  -1/Lg2  000000000000000;... 

0  0  -Rg3/Lg3  0  -1/Lg3  000000000000000;... 

000  -Rg4/Lg4  -1/Lg4  000000000000000;... 

1/Cbus  1/Cbus  1/Cbus  1/Cbus  0  -1/Cbus  -1/Cbus  -1/Cbus  000000000 


0  0  0;  .  .  . 


0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 


1/Lzl  -Rzl/Lzl  0  0  -1/Lzl  00000000000;... 

1/Lz2  0  -Rz2/Lz2  0  0  -1/Lz2  0000000000;... 

1/Lz3  0  0  -Rz3/Lz3  0  0  -1/Lz3  000000000;... 

0  1/Cdl  0  0  -1/Cdl/Rdl  0  0  1/Cdl/Rdl  0  0  -dl/Cdl  0  0  0  0  0;... 

0  0  1/ Cd2  0  0  -1/Cd2/Rd2  0  0  1/Cd2/Rd2  0  0  -d2/Cd2  0  0  0  0;... 

000  1/Cd3  0  0  -1/Cd3/Rd3  0  0  1/Cd3/Rd3  0  0  -d3/Cd3  0  0  0;... 

0000  1/Cdl/Rdl  0  0  -1/Cdl/Rdl  00000000;... 

00000  1/ Cd2/Rd2  0  0  -1/Cd2/Rd2  0000000;... 
000000  1/Cd3/Rd3  0  0  -1/Cd3/Rd3  000000;... 

0000  dl/Lbl  00000000  -1/Lbl  0  0; . . . 

00000  d2/Lb2  00000000  -1/Lb2  0; . . . 

0  0  0  0  0  0  d3/Lb3  00000000  -1/Lb3;  . .  . 
0000000000  1/Cbl  0  0  1/ (Cbl*Rnll )  0  0;... 
00000000000  1/Cb2  0  0  l/(Cb2*Rnl2)  0;... 
000000000000  1/Cb3  0  0  1/ (Cb3*Rnl3) ] ; 


B  =  [Bll/Lgl  0000000000000000000;... 
0  B22/Lg2  000000000000000000;... 

0  0  B33/Lg3  00000000000000000;... 

000  B44/Lg4  0000000000000000;... 

00000000000000000000;... 

00000000000000000000;... 
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0]  ; 

Q 

= 

[Qi  (1) 

0 

0 

0 

0 

0 

0  C 

0 

0 

0 

0 

0 

0 

0 

0  0 

0 
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0000000000000000 

0000000000000000 

0000000000000000 


0  Ri ( 4 )  0  0;..  .%Isl  0.1  0.25 
0  0  Ri (5)  0; . . .%Is2  0.3  0.50 
000  Pen] ;  %Input  penalty  matrix 


try 

[XX,  LL,  GG]  =  care (A,  B,  Q,  R_) ;  %Riccati  Solver 
catch 

GG  =  zeros ( size (A) ) ; 
end 


U  =  -GG*X;  %Input  vector 
end 


8.  Harmonic  Side-Lobe  Detector 

%%%Sidelobe  detection  routine%%% 
function  [SLL]  =  sldetect (Q) 

%  %Get  the  maginitude  FFT  of  the  desired  signal 
mag  =  abs(fft(Q)); 

%normalize  the  FFT 

norm  =  mag. / (max (mag) )  ; 

db  =  20*logl0 (norm) ; 

if  length (db)  >  100 

SLL  =  max (db (30 : 100) ) ; 

else 

SLL  =  1; 

end 

end 


120 


APPENDIX  B.  SEED  CANDIDATE  CALCULATIONS 


%Circuit  Data  for  case  Study  25 
%Circuit  Parameters 


Vref  = 

^  12e3 ; 

%Main 

Bus  voltage 

Vref  1 

=  le3 ; 

%Zone 

1  voltage 

Vref  2 

=  6e3 ; 

%Zone 

2  voltage 

Vref  3 

=  10e3; 

%Zone 

3  voltage 

Rgl  =  4*62. 5e-3;  %Ohm 

Lgl  =  4*500e-6;  %Henry 

Rg2  =  4*75e-3;%Ohm 

Lg2  =  4*450e-6; %Henry 

Rg3  =  4.1*62.5e-3;  %Ohm 

Lg3  =  3.9*500e-6;  %Henry 

Rg4  =  4 . 2*75e-3; %Ohm 

Lg4  =  3 . 8*450e-6; %Henry 

Cbus  =  500e-5;  %6 . 5e-3; %Farad 

P2  =  [20  30  46]*le6; 

PI  =  [15  5  28] *le6; 

fswitch_gen  =  4.0e3; 

Tgen  =  1/f switch_gen; 
fswitch_ESD  =  4*f switch_gen; 

Tesd  =  1/f switch_ESD; 

Fsw_buck  =  le3; 

%Cable  Variables 

%Based  on  the  Cupelli,  de  Carro,  Monti  2015  paper 

%For  800mm  cross  section  cable,  rated  for  871  Amps,  22.1uOhm/m, 

470nH/m, 

%4100pF/m 

%8000A/871  =>  10  parallel  cables 
rho  =  22.1e-6;  %per-meter  reistance 
ind  =  470e-9;  %per  meter  inductance 
cap  =  4100e-12;  %per  meter  capacitance 

%Zone#l-20MW  Ship's  Service 
N1  =  ceil (P2 (1) /Vref/871) ; 

Len  =  300;  %length  of  cable  in  meters 
Rzl  =  Len*rho/Nl;  %Ohm 
Lzl  =  Len*ind/Nl;  %Henry 
Czl  =  Len*cap*Nl;  %Faraday 

dl  =  Vrefl/Vref; 

Reql  =  Vref 1 A2/P1  (1) ; 

Lbl  =  Reql/ (2*Fsw_buck) * (1-dl) ; 

Cripplel  =  (1-dl) /Fsw_buckA2/ (8*Lbl* . 05) ; 


%Zone#2-30MW  Pulsed  Loads 
N2  =  ceil (P2 (2) /Vref/871) ; 
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Len  =  300;  %length  of  cable  in  meters 
Rz2  =  Len*rho/N2;  %0hm 
Lz2  =  Len*ind/N2;  %Henry 
Cz2  =  Len*cap*N2;  %Faraday 

d2  =  Vref2/Vref; 

Req2  =  Vref 2 A2/P1  (2 )  ; 

Lb2  =  Req2/ (2*Fsw_buck) * (l-d2) ; 

Cripple2  =  (l-d2) /Fsw_buckA2/ (8*Lb2* . 05) ; 


%Zone#3-80MW  Populsion 
N3  =  ceil (P2 (3) /Vref/871) ; 

Len  =  300;  %length  of  cable  in  meters 
Rz3  =  Len*rho/N3;  %Ohm 
Lz3  =  Len*ind/N3;  %Henry 
Cz3  =  Len*cap*N3;  %Faraday 

d3  =  Vref3/Vref; 

Req3  =  Vref3A2/Pl  (3) ; 

Lb3  =  Req3/ (2*Fsw_buck) * (l-d3) ; 

Cripple3  =  (l-d3) /Fsw_buckA2/ (8*Lb3* . 05) ; 


taubus  =  max([Lgl/Rgl  Lg2/Rg2  Lg3/Rg3  Lg4/Rg4]); 
Reqbus  =  Vref A2/sum (P2 ) ; 

Cbus  =  taubus/Reqbus ;  %0.0039;  %Farad 
Cdl  =  P2 ( 1 ) /sum (P2 ) *Cbus ; 

Cd2  =  P2 (2) /sum (P2) *Cbus; 

Cd3  =  P2 (3) /sum (P2) *Cbus; 


Rg  =  ( 1/Rgl  +  1/Rg2  +  1/Rg3  +  l/Rg4)A-l; 
Z1  =  (l/dlA2) *VreflA2/P2 (1)  +  Rzl; 

Z2  =  (l/d2A2) *Vref2A2/P2 (2)  +  Rz2; 

Z3  =  (l/d3A2) *Vref3A2/P2 (3)  +  Rz3; 


Zinl 

Zin2 

Zin3 


(dlA2) * (  Rzl  + 
(d2A2) *  (  Rz2  + 
(d3A2) * (  Rz3  + 


(RgA-l  +  Z2A-1  +  Z3A-1)A-1  ) 
(RgA-l  +  Z1A-1  +  Z3A-1)A-1  ) 
(RgA-l  +  Z1A-1  +  Z2A-1)A-1  ) 


Cl  =  Lbl/Reql/ Zinl; 
C2  =  Lb2/Req2/Zin2 ; 
C3  =  Lb3/Req3/Zin3; 


Cbl  =  max (Cripplel ,  Cl); 
Cb2  =  max (Cripple2 ,  C2); 
Cb3  =  max (Cripple3,  4*C3); 


Cx  =  [Cbus  Cdl  Cd2  Cd3  Cbl  Cb2  Cb3]; 


z  =  10; 

Rdl  =  z ;  Rd2  =  z ;  Rd3  =  z ; 

C  =  [Cbus  Cdl  Cd2  Cd3  Cbl  Cb2  Cb3]; 
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